1/* 2 * Copyright 2010 INRIA Saclay 3 * 4 * Use of this software is governed by the MIT license 5 * 6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France, 7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod, 8 * 91893 Orsay, France 9 */ 10 11#include <stdlib.h> 12#include <isl_ctx_private.h> 13#include <isl_map_private.h> 14#include <isl_factorization.h> 15#include <isl_lp_private.h> 16#include <isl_seq.h> 17#include <isl_union_map_private.h> 18#include <isl_constraint_private.h> 19#include <isl_polynomial_private.h> 20#include <isl_point_private.h> 21#include <isl_space_private.h> 22#include <isl_mat_private.h> 23#include <isl_vec_private.h> 24#include <isl_range.h> 25#include <isl_local.h> 26#include <isl_local_space_private.h> 27#include <isl_aff_private.h> 28#include <isl_val_private.h> 29#include <isl_config.h> 30 31#undef EL_BASE 32#define EL_BASE qpolynomial 33 34#include <isl_list_templ.c> 35 36#undef EL_BASE 37#define EL_BASE pw_qpolynomial 38 39#include <isl_list_templ.c> 40 41static unsigned pos(__isl_keep isl_space *space, enum isl_dim_type type) 42{ 43 switch (type) { 44 case isl_dim_param: return 0; 45 case isl_dim_in: return space->nparam; 46 case isl_dim_out: return space->nparam + space->n_in; 47 default: return 0; 48 } 49} 50 51isl_bool isl_poly_is_cst(__isl_keep isl_poly *poly) 52{ 53 if (!poly) 54 return isl_bool_error; 55 56 return isl_bool_ok(poly->var < 0); 57} 58 59__isl_keep isl_poly_cst *isl_poly_as_cst(__isl_keep isl_poly *poly) 60{ 61 if (!poly) 62 return NULL; 63 64 isl_assert(poly->ctx, poly->var < 0, return NULL); 65 66 return (isl_poly_cst *) poly; 67} 68 69__isl_keep isl_poly_rec *isl_poly_as_rec(__isl_keep isl_poly *poly) 70{ 71 if (!poly) 72 return NULL; 73 74 isl_assert(poly->ctx, poly->var >= 0, return NULL); 75 76 return (isl_poly_rec *) poly; 77} 78 79/* Compare two polynomials. 80 * 81 * Return -1 if "poly1" is "smaller" than "poly2", 1 if "poly1" is "greater" 82 * than "poly2" and 0 if they are equal. 83 */ 84static int isl_poly_plain_cmp(__isl_keep isl_poly *poly1, 85 __isl_keep isl_poly *poly2) 86{ 87 int i; 88 isl_bool is_cst1; 89 isl_poly_rec *rec1, *rec2; 90 91 if (poly1 == poly2) 92 return 0; 93 is_cst1 = isl_poly_is_cst(poly1); 94 if (is_cst1 < 0) 95 return -1; 96 if (!poly2) 97 return 1; 98 if (poly1->var != poly2->var) 99 return poly1->var - poly2->var; 100 101 if (is_cst1) { 102 isl_poly_cst *cst1, *cst2; 103 int cmp; 104 105 cst1 = isl_poly_as_cst(poly1); 106 cst2 = isl_poly_as_cst(poly2); 107 if (!cst1 || !cst2) 108 return 0; 109 cmp = isl_int_cmp(cst1->n, cst2->n); 110 if (cmp != 0) 111 return cmp; 112 return isl_int_cmp(cst1->d, cst2->d); 113 } 114 115 rec1 = isl_poly_as_rec(poly1); 116 rec2 = isl_poly_as_rec(poly2); 117 if (!rec1 || !rec2) 118 return 0; 119 120 if (rec1->n != rec2->n) 121 return rec1->n - rec2->n; 122 123 for (i = 0; i < rec1->n; ++i) { 124 int cmp = isl_poly_plain_cmp(rec1->p[i], rec2->p[i]); 125 if (cmp != 0) 126 return cmp; 127 } 128 129 return 0; 130} 131 132isl_bool isl_poly_is_equal(__isl_keep isl_poly *poly1, 133 __isl_keep isl_poly *poly2) 134{ 135 int i; 136 isl_bool is_cst1; 137 isl_poly_rec *rec1, *rec2; 138 139 is_cst1 = isl_poly_is_cst(poly1); 140 if (is_cst1 < 0 || !poly2) 141 return isl_bool_error; 142 if (poly1 == poly2) 143 return isl_bool_true; 144 if (poly1->var != poly2->var) 145 return isl_bool_false; 146 if (is_cst1) { 147 isl_poly_cst *cst1, *cst2; 148 int r; 149 cst1 = isl_poly_as_cst(poly1); 150 cst2 = isl_poly_as_cst(poly2); 151 if (!cst1 || !cst2) 152 return isl_bool_error; 153 r = isl_int_eq(cst1->n, cst2->n) && 154 isl_int_eq(cst1->d, cst2->d); 155 return isl_bool_ok(r); 156 } 157 158 rec1 = isl_poly_as_rec(poly1); 159 rec2 = isl_poly_as_rec(poly2); 160 if (!rec1 || !rec2) 161 return isl_bool_error; 162 163 if (rec1->n != rec2->n) 164 return isl_bool_false; 165 166 for (i = 0; i < rec1->n; ++i) { 167 isl_bool eq = isl_poly_is_equal(rec1->p[i], rec2->p[i]); 168 if (eq < 0 || !eq) 169 return eq; 170 } 171 172 return isl_bool_true; 173} 174 175isl_bool isl_poly_is_zero(__isl_keep isl_poly *poly) 176{ 177 isl_bool is_cst; 178 isl_poly_cst *cst; 179 180 is_cst = isl_poly_is_cst(poly); 181 if (is_cst < 0 || !is_cst) 182 return is_cst; 183 184 cst = isl_poly_as_cst(poly); 185 if (!cst) 186 return isl_bool_error; 187 188 return isl_bool_ok(isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d)); 189} 190 191int isl_poly_sgn(__isl_keep isl_poly *poly) 192{ 193 isl_bool is_cst; 194 isl_poly_cst *cst; 195 196 is_cst = isl_poly_is_cst(poly); 197 if (is_cst < 0 || !is_cst) 198 return 0; 199 200 cst = isl_poly_as_cst(poly); 201 if (!cst) 202 return 0; 203 204 return isl_int_sgn(cst->n); 205} 206 207isl_bool isl_poly_is_nan(__isl_keep isl_poly *poly) 208{ 209 isl_bool is_cst; 210 isl_poly_cst *cst; 211 212 is_cst = isl_poly_is_cst(poly); 213 if (is_cst < 0 || !is_cst) 214 return is_cst; 215 216 cst = isl_poly_as_cst(poly); 217 if (!cst) 218 return isl_bool_error; 219 220 return isl_bool_ok(isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d)); 221} 222 223isl_bool isl_poly_is_infty(__isl_keep isl_poly *poly) 224{ 225 isl_bool is_cst; 226 isl_poly_cst *cst; 227 228 is_cst = isl_poly_is_cst(poly); 229 if (is_cst < 0 || !is_cst) 230 return is_cst; 231 232 cst = isl_poly_as_cst(poly); 233 if (!cst) 234 return isl_bool_error; 235 236 return isl_bool_ok(isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d)); 237} 238 239isl_bool isl_poly_is_neginfty(__isl_keep isl_poly *poly) 240{ 241 isl_bool is_cst; 242 isl_poly_cst *cst; 243 244 is_cst = isl_poly_is_cst(poly); 245 if (is_cst < 0 || !is_cst) 246 return is_cst; 247 248 cst = isl_poly_as_cst(poly); 249 if (!cst) 250 return isl_bool_error; 251 252 return isl_bool_ok(isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d)); 253} 254 255isl_bool isl_poly_is_one(__isl_keep isl_poly *poly) 256{ 257 isl_bool is_cst; 258 isl_poly_cst *cst; 259 int r; 260 261 is_cst = isl_poly_is_cst(poly); 262 if (is_cst < 0 || !is_cst) 263 return is_cst; 264 265 cst = isl_poly_as_cst(poly); 266 if (!cst) 267 return isl_bool_error; 268 269 r = isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d); 270 return isl_bool_ok(r); 271} 272 273isl_bool isl_poly_is_negone(__isl_keep isl_poly *poly) 274{ 275 isl_bool is_cst; 276 isl_poly_cst *cst; 277 278 is_cst = isl_poly_is_cst(poly); 279 if (is_cst < 0 || !is_cst) 280 return is_cst; 281 282 cst = isl_poly_as_cst(poly); 283 if (!cst) 284 return isl_bool_error; 285 286 return isl_bool_ok(isl_int_is_negone(cst->n) && isl_int_is_one(cst->d)); 287} 288 289__isl_give isl_poly_cst *isl_poly_cst_alloc(isl_ctx *ctx) 290{ 291 isl_poly_cst *cst; 292 293 cst = isl_alloc_type(ctx, struct isl_poly_cst); 294 if (!cst) 295 return NULL; 296 297 cst->poly.ref = 1; 298 cst->poly.ctx = ctx; 299 isl_ctx_ref(ctx); 300 cst->poly.var = -1; 301 302 isl_int_init(cst->n); 303 isl_int_init(cst->d); 304 305 return cst; 306} 307 308__isl_give isl_poly *isl_poly_zero(isl_ctx *ctx) 309{ 310 isl_poly_cst *cst; 311 312 cst = isl_poly_cst_alloc(ctx); 313 if (!cst) 314 return NULL; 315 316 isl_int_set_si(cst->n, 0); 317 isl_int_set_si(cst->d, 1); 318 319 return &cst->poly; 320} 321 322__isl_give isl_poly *isl_poly_one(isl_ctx *ctx) 323{ 324 isl_poly_cst *cst; 325 326 cst = isl_poly_cst_alloc(ctx); 327 if (!cst) 328 return NULL; 329 330 isl_int_set_si(cst->n, 1); 331 isl_int_set_si(cst->d, 1); 332 333 return &cst->poly; 334} 335 336__isl_give isl_poly *isl_poly_infty(isl_ctx *ctx) 337{ 338 isl_poly_cst *cst; 339 340 cst = isl_poly_cst_alloc(ctx); 341 if (!cst) 342 return NULL; 343 344 isl_int_set_si(cst->n, 1); 345 isl_int_set_si(cst->d, 0); 346 347 return &cst->poly; 348} 349 350__isl_give isl_poly *isl_poly_neginfty(isl_ctx *ctx) 351{ 352 isl_poly_cst *cst; 353 354 cst = isl_poly_cst_alloc(ctx); 355 if (!cst) 356 return NULL; 357 358 isl_int_set_si(cst->n, -1); 359 isl_int_set_si(cst->d, 0); 360 361 return &cst->poly; 362} 363 364__isl_give isl_poly *isl_poly_nan(isl_ctx *ctx) 365{ 366 isl_poly_cst *cst; 367 368 cst = isl_poly_cst_alloc(ctx); 369 if (!cst) 370 return NULL; 371 372 isl_int_set_si(cst->n, 0); 373 isl_int_set_si(cst->d, 0); 374 375 return &cst->poly; 376} 377 378__isl_give isl_poly *isl_poly_rat_cst(isl_ctx *ctx, isl_int n, isl_int d) 379{ 380 isl_poly_cst *cst; 381 382 cst = isl_poly_cst_alloc(ctx); 383 if (!cst) 384 return NULL; 385 386 isl_int_set(cst->n, n); 387 isl_int_set(cst->d, d); 388 389 return &cst->poly; 390} 391 392__isl_give isl_poly_rec *isl_poly_alloc_rec(isl_ctx *ctx, int var, int size) 393{ 394 isl_poly_rec *rec; 395 396 isl_assert(ctx, var >= 0, return NULL); 397 isl_assert(ctx, size >= 0, return NULL); 398 rec = isl_calloc(ctx, struct isl_poly_rec, 399 sizeof(struct isl_poly_rec) + 400 size * sizeof(struct isl_poly *)); 401 if (!rec) 402 return NULL; 403 404 rec->poly.ref = 1; 405 rec->poly.ctx = ctx; 406 isl_ctx_ref(ctx); 407 rec->poly.var = var; 408 409 rec->n = 0; 410 rec->size = size; 411 412 return rec; 413} 414 415/* Return the domain space of "qp". 416 * This may be either a copy or the space itself 417 * if there is only one reference to "qp". 418 * This allows the space to be modified inplace 419 * if both the quasi-polynomial and its domain space 420 * have only a single reference. 421 * The caller is not allowed to modify "qp" between this call and 422 * a subsequent call to isl_qpolynomial_restore_domain_space. 423 * The only exception is that isl_qpolynomial_free can be called instead. 424 */ 425static __isl_give isl_space *isl_qpolynomial_take_domain_space( 426 __isl_keep isl_qpolynomial *qp) 427{ 428 isl_space *space; 429 430 if (!qp) 431 return NULL; 432 if (qp->ref != 1) 433 return isl_qpolynomial_get_domain_space(qp); 434 space = qp->dim; 435 qp->dim = NULL; 436 return space; 437} 438 439/* Set the domain space of "qp" to "space", 440 * where the domain space of "qp" may be missing 441 * due to a preceding call to isl_qpolynomial_take_domain_space. 442 * However, in this case, "qp" only has a single reference and 443 * then the call to isl_qpolynomial_cow has no effect. 444 */ 445static __isl_give isl_qpolynomial *isl_qpolynomial_restore_domain_space( 446 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space) 447{ 448 if (!qp || !space) 449 goto error; 450 451 if (qp->dim == space) { 452 isl_space_free(space); 453 return qp; 454 } 455 456 qp = isl_qpolynomial_cow(qp); 457 if (!qp) 458 goto error; 459 isl_space_free(qp->dim); 460 qp->dim = space; 461 462 return qp; 463error: 464 isl_qpolynomial_free(qp); 465 isl_space_free(space); 466 return NULL; 467} 468 469__isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space( 470 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space) 471{ 472 return isl_qpolynomial_restore_domain_space(qp, space); 473} 474 475/* Reset the space of "qp". This function is called from isl_pw_templ.c 476 * and doesn't know if the space of an element object is represented 477 * directly or through its domain. It therefore passes along both. 478 */ 479__isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain( 480 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space, 481 __isl_take isl_space *domain) 482{ 483 isl_space_free(space); 484 return isl_qpolynomial_reset_domain_space(qp, domain); 485} 486 487isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp) 488{ 489 return qp ? qp->dim->ctx : NULL; 490} 491 492/* Return the domain space of "qp". 493 */ 494static __isl_keep isl_space *isl_qpolynomial_peek_domain_space( 495 __isl_keep isl_qpolynomial *qp) 496{ 497 return qp ? qp->dim : NULL; 498} 499 500/* Return a copy of the domain space of "qp". 501 */ 502__isl_give isl_space *isl_qpolynomial_get_domain_space( 503 __isl_keep isl_qpolynomial *qp) 504{ 505 return isl_space_copy(isl_qpolynomial_peek_domain_space(qp)); 506} 507 508#undef TYPE 509#define TYPE isl_qpolynomial 510#undef PEEK_SPACE 511#define PEEK_SPACE peek_domain_space 512 513static 514#include "isl_type_has_equal_space_bin_templ.c" 515static 516#include "isl_type_check_equal_space_templ.c" 517 518#undef PEEK_SPACE 519 520/* Return a copy of the local variables of "qp". 521 */ 522__isl_keep isl_local *isl_qpolynomial_get_local( 523 __isl_keep isl_qpolynomial *qp) 524{ 525 return qp ? isl_local_copy(qp->div) : NULL; 526} 527 528/* Return the local variables of "qp". 529 * This may be either a copy or the local variables themselves 530 * if there is only one reference to "qp". 531 * This allows the local variables to be modified in-place 532 * if both the quasi-polynomial and its local variables 533 * have only a single reference. 534 * The caller is not allowed to modify "qp" between this call and 535 * the subsequent call to isl_qpolynomial_restore_local. 536 * The only exception is that isl_qpolynomial_free can be called instead. 537 */ 538static __isl_give isl_local *isl_qpolynomial_take_local( 539 __isl_keep isl_qpolynomial *qp) 540{ 541 isl_local *local; 542 543 if (!qp) 544 return NULL; 545 if (qp->ref != 1) 546 return isl_qpolynomial_get_local(qp); 547 local = qp->div; 548 qp->div = NULL; 549 return local; 550} 551 552/* Set the local variables of "qp" to "local", 553 * where the local variables of "qp" may be missing 554 * due to a preceding call to isl_qpolynomial_take_local. 555 * However, in this case, "qp" only has a single reference and 556 * then the call to isl_qpolynomial_cow has no effect. 557 */ 558static __isl_give isl_qpolynomial *isl_qpolynomial_restore_local( 559 __isl_keep isl_qpolynomial *qp, __isl_take isl_local *local) 560{ 561 if (!qp || !local) 562 goto error; 563 564 if (qp->div == local) { 565 isl_local_free(local); 566 return qp; 567 } 568 569 qp = isl_qpolynomial_cow(qp); 570 if (!qp) 571 goto error; 572 isl_local_free(qp->div); 573 qp->div = local; 574 575 return qp; 576error: 577 isl_qpolynomial_free(qp); 578 isl_local_free(local); 579 return NULL; 580} 581 582/* Return a copy of the local space on which "qp" is defined. 583 */ 584static __isl_give isl_local_space *isl_qpolynomial_get_domain_local_space( 585 __isl_keep isl_qpolynomial *qp) 586{ 587 isl_space *space; 588 isl_local *local; 589 590 if (!qp) 591 return NULL; 592 593 space = isl_qpolynomial_get_domain_space(qp); 594 local = isl_qpolynomial_get_local(qp); 595 return isl_local_space_alloc_div(space, local); 596} 597 598__isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp) 599{ 600 isl_space *space; 601 if (!qp) 602 return NULL; 603 space = isl_space_copy(qp->dim); 604 space = isl_space_from_domain(space); 605 space = isl_space_add_dims(space, isl_dim_out, 1); 606 return space; 607} 608 609/* Return the number of variables of the given type in the domain of "qp". 610 */ 611isl_size isl_qpolynomial_domain_dim(__isl_keep isl_qpolynomial *qp, 612 enum isl_dim_type type) 613{ 614 isl_space *space; 615 isl_size dim; 616 617 space = isl_qpolynomial_peek_domain_space(qp); 618 619 if (!space) 620 return isl_size_error; 621 if (type == isl_dim_div) 622 return qp->div->n_row; 623 dim = isl_space_dim(space, type); 624 if (dim < 0) 625 return isl_size_error; 626 if (type == isl_dim_all) { 627 isl_size n_div; 628 629 n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div); 630 if (n_div < 0) 631 return isl_size_error; 632 dim += n_div; 633 } 634 return dim; 635} 636 637/* Given the type of a dimension of an isl_qpolynomial, 638 * return the type of the corresponding dimension in its domain. 639 * This function is only called for "type" equal to isl_dim_in or 640 * isl_dim_param. 641 */ 642static enum isl_dim_type domain_type(enum isl_dim_type type) 643{ 644 return type == isl_dim_in ? isl_dim_set : type; 645} 646 647/* Externally, an isl_qpolynomial has a map space, but internally, the 648 * ls field corresponds to the domain of that space. 649 */ 650isl_size isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp, 651 enum isl_dim_type type) 652{ 653 if (!qp) 654 return isl_size_error; 655 if (type == isl_dim_out) 656 return 1; 657 type = domain_type(type); 658 return isl_qpolynomial_domain_dim(qp, type); 659} 660 661/* Return the offset of the first variable of type "type" within 662 * the variables of the domain of "qp". 663 */ 664static isl_size isl_qpolynomial_domain_var_offset( 665 __isl_keep isl_qpolynomial *qp, enum isl_dim_type type) 666{ 667 isl_space *space; 668 669 space = isl_qpolynomial_peek_domain_space(qp); 670 671 switch (type) { 672 case isl_dim_param: 673 case isl_dim_set: return isl_space_offset(space, type); 674 case isl_dim_div: return isl_space_dim(space, isl_dim_all); 675 case isl_dim_cst: 676 default: 677 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid, 678 "invalid dimension type", return isl_size_error); 679 } 680} 681 682/* Return the offset of the first coefficient of type "type" in 683 * the domain of "qp". 684 */ 685unsigned isl_qpolynomial_domain_offset(__isl_keep isl_qpolynomial *qp, 686 enum isl_dim_type type) 687{ 688 switch (type) { 689 case isl_dim_cst: 690 return 0; 691 case isl_dim_param: 692 case isl_dim_set: 693 case isl_dim_div: 694 return 1 + isl_qpolynomial_domain_var_offset(qp, type); 695 default: 696 return 0; 697 } 698} 699 700isl_bool isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp) 701{ 702 return qp ? isl_poly_is_zero(qp->poly) : isl_bool_error; 703} 704 705isl_bool isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp) 706{ 707 return qp ? isl_poly_is_one(qp->poly) : isl_bool_error; 708} 709 710isl_bool isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp) 711{ 712 return qp ? isl_poly_is_nan(qp->poly) : isl_bool_error; 713} 714 715isl_bool isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp) 716{ 717 return qp ? isl_poly_is_infty(qp->poly) : isl_bool_error; 718} 719 720isl_bool isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp) 721{ 722 return qp ? isl_poly_is_neginfty(qp->poly) : isl_bool_error; 723} 724 725int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp) 726{ 727 return qp ? isl_poly_sgn(qp->poly) : 0; 728} 729 730static void poly_free_cst(__isl_take isl_poly_cst *cst) 731{ 732 isl_int_clear(cst->n); 733 isl_int_clear(cst->d); 734} 735 736static void poly_free_rec(__isl_take isl_poly_rec *rec) 737{ 738 int i; 739 740 for (i = 0; i < rec->n; ++i) 741 isl_poly_free(rec->p[i]); 742} 743 744__isl_give isl_poly *isl_poly_copy(__isl_keep isl_poly *poly) 745{ 746 if (!poly) 747 return NULL; 748 749 poly->ref++; 750 return poly; 751} 752 753__isl_give isl_poly *isl_poly_dup_cst(__isl_keep isl_poly *poly) 754{ 755 isl_poly_cst *cst; 756 isl_poly_cst *dup; 757 758 cst = isl_poly_as_cst(poly); 759 if (!cst) 760 return NULL; 761 762 dup = isl_poly_as_cst(isl_poly_zero(poly->ctx)); 763 if (!dup) 764 return NULL; 765 isl_int_set(dup->n, cst->n); 766 isl_int_set(dup->d, cst->d); 767 768 return &dup->poly; 769} 770 771__isl_give isl_poly *isl_poly_dup_rec(__isl_keep isl_poly *poly) 772{ 773 int i; 774 isl_poly_rec *rec; 775 isl_poly_rec *dup; 776 777 rec = isl_poly_as_rec(poly); 778 if (!rec) 779 return NULL; 780 781 dup = isl_poly_alloc_rec(poly->ctx, poly->var, rec->n); 782 if (!dup) 783 return NULL; 784 785 for (i = 0; i < rec->n; ++i) { 786 dup->p[i] = isl_poly_copy(rec->p[i]); 787 if (!dup->p[i]) 788 goto error; 789 dup->n++; 790 } 791 792 return &dup->poly; 793error: 794 isl_poly_free(&dup->poly); 795 return NULL; 796} 797 798__isl_give isl_poly *isl_poly_dup(__isl_keep isl_poly *poly) 799{ 800 isl_bool is_cst; 801 802 is_cst = isl_poly_is_cst(poly); 803 if (is_cst < 0) 804 return NULL; 805 if (is_cst) 806 return isl_poly_dup_cst(poly); 807 else 808 return isl_poly_dup_rec(poly); 809} 810 811__isl_give isl_poly *isl_poly_cow(__isl_take isl_poly *poly) 812{ 813 if (!poly) 814 return NULL; 815 816 if (poly->ref == 1) 817 return poly; 818 poly->ref--; 819 return isl_poly_dup(poly); 820} 821 822__isl_null isl_poly *isl_poly_free(__isl_take isl_poly *poly) 823{ 824 if (!poly) 825 return NULL; 826 827 if (--poly->ref > 0) 828 return NULL; 829 830 if (poly->var < 0) 831 poly_free_cst((isl_poly_cst *) poly); 832 else 833 poly_free_rec((isl_poly_rec *) poly); 834 835 isl_ctx_deref(poly->ctx); 836 free(poly); 837 return NULL; 838} 839 840static void isl_poly_cst_reduce(__isl_keep isl_poly_cst *cst) 841{ 842 isl_int gcd; 843 844 isl_int_init(gcd); 845 isl_int_gcd(gcd, cst->n, cst->d); 846 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) { 847 isl_int_divexact(cst->n, cst->n, gcd); 848 isl_int_divexact(cst->d, cst->d, gcd); 849 } 850 isl_int_clear(gcd); 851} 852 853__isl_give isl_poly *isl_poly_sum_cst(__isl_take isl_poly *poly1, 854 __isl_take isl_poly *poly2) 855{ 856 isl_poly_cst *cst1; 857 isl_poly_cst *cst2; 858 859 poly1 = isl_poly_cow(poly1); 860 if (!poly1 || !poly2) 861 goto error; 862 863 cst1 = isl_poly_as_cst(poly1); 864 cst2 = isl_poly_as_cst(poly2); 865 866 if (isl_int_eq(cst1->d, cst2->d)) 867 isl_int_add(cst1->n, cst1->n, cst2->n); 868 else { 869 isl_int_mul(cst1->n, cst1->n, cst2->d); 870 isl_int_addmul(cst1->n, cst2->n, cst1->d); 871 isl_int_mul(cst1->d, cst1->d, cst2->d); 872 } 873 874 isl_poly_cst_reduce(cst1); 875 876 isl_poly_free(poly2); 877 return poly1; 878error: 879 isl_poly_free(poly1); 880 isl_poly_free(poly2); 881 return NULL; 882} 883 884static __isl_give isl_poly *replace_by_zero(__isl_take isl_poly *poly) 885{ 886 struct isl_ctx *ctx; 887 888 if (!poly) 889 return NULL; 890 ctx = poly->ctx; 891 isl_poly_free(poly); 892 return isl_poly_zero(ctx); 893} 894 895static __isl_give isl_poly *replace_by_constant_term(__isl_take isl_poly *poly) 896{ 897 isl_poly_rec *rec; 898 isl_poly *cst; 899 900 if (!poly) 901 return NULL; 902 903 rec = isl_poly_as_rec(poly); 904 if (!rec) 905 goto error; 906 cst = isl_poly_copy(rec->p[0]); 907 isl_poly_free(poly); 908 return cst; 909error: 910 isl_poly_free(poly); 911 return NULL; 912} 913 914__isl_give isl_poly *isl_poly_sum(__isl_take isl_poly *poly1, 915 __isl_take isl_poly *poly2) 916{ 917 int i; 918 isl_bool is_zero, is_nan, is_cst; 919 isl_poly_rec *rec1, *rec2; 920 921 if (!poly1 || !poly2) 922 goto error; 923 924 is_nan = isl_poly_is_nan(poly1); 925 if (is_nan < 0) 926 goto error; 927 if (is_nan) { 928 isl_poly_free(poly2); 929 return poly1; 930 } 931 932 is_nan = isl_poly_is_nan(poly2); 933 if (is_nan < 0) 934 goto error; 935 if (is_nan) { 936 isl_poly_free(poly1); 937 return poly2; 938 } 939 940 is_zero = isl_poly_is_zero(poly1); 941 if (is_zero < 0) 942 goto error; 943 if (is_zero) { 944 isl_poly_free(poly1); 945 return poly2; 946 } 947 948 is_zero = isl_poly_is_zero(poly2); 949 if (is_zero < 0) 950 goto error; 951 if (is_zero) { 952 isl_poly_free(poly2); 953 return poly1; 954 } 955 956 if (poly1->var < poly2->var) 957 return isl_poly_sum(poly2, poly1); 958 959 if (poly2->var < poly1->var) { 960 isl_poly_rec *rec; 961 isl_bool is_infty; 962 963 is_infty = isl_poly_is_infty(poly2); 964 if (is_infty >= 0 && !is_infty) 965 is_infty = isl_poly_is_neginfty(poly2); 966 if (is_infty < 0) 967 goto error; 968 if (is_infty) { 969 isl_poly_free(poly1); 970 return poly2; 971 } 972 poly1 = isl_poly_cow(poly1); 973 rec = isl_poly_as_rec(poly1); 974 if (!rec) 975 goto error; 976 rec->p[0] = isl_poly_sum(rec->p[0], poly2); 977 if (rec->n == 1) 978 poly1 = replace_by_constant_term(poly1); 979 return poly1; 980 } 981 982 is_cst = isl_poly_is_cst(poly1); 983 if (is_cst < 0) 984 goto error; 985 if (is_cst) 986 return isl_poly_sum_cst(poly1, poly2); 987 988 rec1 = isl_poly_as_rec(poly1); 989 rec2 = isl_poly_as_rec(poly2); 990 if (!rec1 || !rec2) 991 goto error; 992 993 if (rec1->n < rec2->n) 994 return isl_poly_sum(poly2, poly1); 995 996 poly1 = isl_poly_cow(poly1); 997 rec1 = isl_poly_as_rec(poly1); 998 if (!rec1) 999 goto error; 1000 1001 for (i = rec2->n - 1; i >= 0; --i) { 1002 isl_bool is_zero; 1003 1004 rec1->p[i] = isl_poly_sum(rec1->p[i], 1005 isl_poly_copy(rec2->p[i])); 1006 if (!rec1->p[i]) 1007 goto error; 1008 if (i != rec1->n - 1) 1009 continue; 1010 is_zero = isl_poly_is_zero(rec1->p[i]); 1011 if (is_zero < 0) 1012 goto error; 1013 if (is_zero) { 1014 isl_poly_free(rec1->p[i]); 1015 rec1->n--; 1016 } 1017 } 1018 1019 if (rec1->n == 0) 1020 poly1 = replace_by_zero(poly1); 1021 else if (rec1->n == 1) 1022 poly1 = replace_by_constant_term(poly1); 1023 1024 isl_poly_free(poly2); 1025 1026 return poly1; 1027error: 1028 isl_poly_free(poly1); 1029 isl_poly_free(poly2); 1030 return NULL; 1031} 1032 1033__isl_give isl_poly *isl_poly_cst_add_isl_int(__isl_take isl_poly *poly, 1034 isl_int v) 1035{ 1036 isl_poly_cst *cst; 1037 1038 poly = isl_poly_cow(poly); 1039 if (!poly) 1040 return NULL; 1041 1042 cst = isl_poly_as_cst(poly); 1043 1044 isl_int_addmul(cst->n, cst->d, v); 1045 1046 return poly; 1047} 1048 1049__isl_give isl_poly *isl_poly_add_isl_int(__isl_take isl_poly *poly, isl_int v) 1050{ 1051 isl_bool is_cst; 1052 isl_poly_rec *rec; 1053 1054 is_cst = isl_poly_is_cst(poly); 1055 if (is_cst < 0) 1056 return isl_poly_free(poly); 1057 if (is_cst) 1058 return isl_poly_cst_add_isl_int(poly, v); 1059 1060 poly = isl_poly_cow(poly); 1061 rec = isl_poly_as_rec(poly); 1062 if (!rec) 1063 goto error; 1064 1065 rec->p[0] = isl_poly_add_isl_int(rec->p[0], v); 1066 if (!rec->p[0]) 1067 goto error; 1068 1069 return poly; 1070error: 1071 isl_poly_free(poly); 1072 return NULL; 1073} 1074 1075__isl_give isl_poly *isl_poly_cst_mul_isl_int(__isl_take isl_poly *poly, 1076 isl_int v) 1077{ 1078 isl_bool is_zero; 1079 isl_poly_cst *cst; 1080 1081 is_zero = isl_poly_is_zero(poly); 1082 if (is_zero < 0) 1083 return isl_poly_free(poly); 1084 if (is_zero) 1085 return poly; 1086 1087 poly = isl_poly_cow(poly); 1088 if (!poly) 1089 return NULL; 1090 1091 cst = isl_poly_as_cst(poly); 1092 1093 isl_int_mul(cst->n, cst->n, v); 1094 1095 return poly; 1096} 1097 1098__isl_give isl_poly *isl_poly_mul_isl_int(__isl_take isl_poly *poly, isl_int v) 1099{ 1100 int i; 1101 isl_bool is_cst; 1102 isl_poly_rec *rec; 1103 1104 is_cst = isl_poly_is_cst(poly); 1105 if (is_cst < 0) 1106 return isl_poly_free(poly); 1107 if (is_cst) 1108 return isl_poly_cst_mul_isl_int(poly, v); 1109 1110 poly = isl_poly_cow(poly); 1111 rec = isl_poly_as_rec(poly); 1112 if (!rec) 1113 goto error; 1114 1115 for (i = 0; i < rec->n; ++i) { 1116 rec->p[i] = isl_poly_mul_isl_int(rec->p[i], v); 1117 if (!rec->p[i]) 1118 goto error; 1119 } 1120 1121 return poly; 1122error: 1123 isl_poly_free(poly); 1124 return NULL; 1125} 1126 1127/* Multiply the constant polynomial "poly" by "v". 1128 */ 1129static __isl_give isl_poly *isl_poly_cst_scale_val(__isl_take isl_poly *poly, 1130 __isl_keep isl_val *v) 1131{ 1132 isl_bool is_zero; 1133 isl_poly_cst *cst; 1134 1135 is_zero = isl_poly_is_zero(poly); 1136 if (is_zero < 0) 1137 return isl_poly_free(poly); 1138 if (is_zero) 1139 return poly; 1140 1141 poly = isl_poly_cow(poly); 1142 if (!poly) 1143 return NULL; 1144 1145 cst = isl_poly_as_cst(poly); 1146 1147 isl_int_mul(cst->n, cst->n, v->n); 1148 isl_int_mul(cst->d, cst->d, v->d); 1149 isl_poly_cst_reduce(cst); 1150 1151 return poly; 1152} 1153 1154/* Multiply the polynomial "poly" by "v". 1155 */ 1156static __isl_give isl_poly *isl_poly_scale_val(__isl_take isl_poly *poly, 1157 __isl_keep isl_val *v) 1158{ 1159 int i; 1160 isl_bool is_cst; 1161 isl_poly_rec *rec; 1162 1163 is_cst = isl_poly_is_cst(poly); 1164 if (is_cst < 0) 1165 return isl_poly_free(poly); 1166 if (is_cst) 1167 return isl_poly_cst_scale_val(poly, v); 1168 1169 poly = isl_poly_cow(poly); 1170 rec = isl_poly_as_rec(poly); 1171 if (!rec) 1172 goto error; 1173 1174 for (i = 0; i < rec->n; ++i) { 1175 rec->p[i] = isl_poly_scale_val(rec->p[i], v); 1176 if (!rec->p[i]) 1177 goto error; 1178 } 1179 1180 return poly; 1181error: 1182 isl_poly_free(poly); 1183 return NULL; 1184} 1185 1186__isl_give isl_poly *isl_poly_mul_cst(__isl_take isl_poly *poly1, 1187 __isl_take isl_poly *poly2) 1188{ 1189 isl_poly_cst *cst1; 1190 isl_poly_cst *cst2; 1191 1192 poly1 = isl_poly_cow(poly1); 1193 if (!poly1 || !poly2) 1194 goto error; 1195 1196 cst1 = isl_poly_as_cst(poly1); 1197 cst2 = isl_poly_as_cst(poly2); 1198 1199 isl_int_mul(cst1->n, cst1->n, cst2->n); 1200 isl_int_mul(cst1->d, cst1->d, cst2->d); 1201 1202 isl_poly_cst_reduce(cst1); 1203 1204 isl_poly_free(poly2); 1205 return poly1; 1206error: 1207 isl_poly_free(poly1); 1208 isl_poly_free(poly2); 1209 return NULL; 1210} 1211 1212__isl_give isl_poly *isl_poly_mul_rec(__isl_take isl_poly *poly1, 1213 __isl_take isl_poly *poly2) 1214{ 1215 isl_poly_rec *rec1; 1216 isl_poly_rec *rec2; 1217 isl_poly_rec *res = NULL; 1218 int i, j; 1219 int size; 1220 1221 rec1 = isl_poly_as_rec(poly1); 1222 rec2 = isl_poly_as_rec(poly2); 1223 if (!rec1 || !rec2) 1224 goto error; 1225 size = rec1->n + rec2->n - 1; 1226 res = isl_poly_alloc_rec(poly1->ctx, poly1->var, size); 1227 if (!res) 1228 goto error; 1229 1230 for (i = 0; i < rec1->n; ++i) { 1231 res->p[i] = isl_poly_mul(isl_poly_copy(rec2->p[0]), 1232 isl_poly_copy(rec1->p[i])); 1233 if (!res->p[i]) 1234 goto error; 1235 res->n++; 1236 } 1237 for (; i < size; ++i) { 1238 res->p[i] = isl_poly_zero(poly1->ctx); 1239 if (!res->p[i]) 1240 goto error; 1241 res->n++; 1242 } 1243 for (i = 0; i < rec1->n; ++i) { 1244 for (j = 1; j < rec2->n; ++j) { 1245 isl_poly *poly; 1246 poly = isl_poly_mul(isl_poly_copy(rec2->p[j]), 1247 isl_poly_copy(rec1->p[i])); 1248 res->p[i + j] = isl_poly_sum(res->p[i + j], poly); 1249 if (!res->p[i + j]) 1250 goto error; 1251 } 1252 } 1253 1254 isl_poly_free(poly1); 1255 isl_poly_free(poly2); 1256 1257 return &res->poly; 1258error: 1259 isl_poly_free(poly1); 1260 isl_poly_free(poly2); 1261 isl_poly_free(&res->poly); 1262 return NULL; 1263} 1264 1265__isl_give isl_poly *isl_poly_mul(__isl_take isl_poly *poly1, 1266 __isl_take isl_poly *poly2) 1267{ 1268 isl_bool is_zero, is_nan, is_one, is_cst; 1269 1270 if (!poly1 || !poly2) 1271 goto error; 1272 1273 is_nan = isl_poly_is_nan(poly1); 1274 if (is_nan < 0) 1275 goto error; 1276 if (is_nan) { 1277 isl_poly_free(poly2); 1278 return poly1; 1279 } 1280 1281 is_nan = isl_poly_is_nan(poly2); 1282 if (is_nan < 0) 1283 goto error; 1284 if (is_nan) { 1285 isl_poly_free(poly1); 1286 return poly2; 1287 } 1288 1289 is_zero = isl_poly_is_zero(poly1); 1290 if (is_zero < 0) 1291 goto error; 1292 if (is_zero) { 1293 isl_poly_free(poly2); 1294 return poly1; 1295 } 1296 1297 is_zero = isl_poly_is_zero(poly2); 1298 if (is_zero < 0) 1299 goto error; 1300 if (is_zero) { 1301 isl_poly_free(poly1); 1302 return poly2; 1303 } 1304 1305 is_one = isl_poly_is_one(poly1); 1306 if (is_one < 0) 1307 goto error; 1308 if (is_one) { 1309 isl_poly_free(poly1); 1310 return poly2; 1311 } 1312 1313 is_one = isl_poly_is_one(poly2); 1314 if (is_one < 0) 1315 goto error; 1316 if (is_one) { 1317 isl_poly_free(poly2); 1318 return poly1; 1319 } 1320 1321 if (poly1->var < poly2->var) 1322 return isl_poly_mul(poly2, poly1); 1323 1324 if (poly2->var < poly1->var) { 1325 int i; 1326 isl_poly_rec *rec; 1327 isl_bool is_infty; 1328 1329 is_infty = isl_poly_is_infty(poly2); 1330 if (is_infty >= 0 && !is_infty) 1331 is_infty = isl_poly_is_neginfty(poly2); 1332 if (is_infty < 0) 1333 goto error; 1334 if (is_infty) { 1335 isl_ctx *ctx = poly1->ctx; 1336 isl_poly_free(poly1); 1337 isl_poly_free(poly2); 1338 return isl_poly_nan(ctx); 1339 } 1340 poly1 = isl_poly_cow(poly1); 1341 rec = isl_poly_as_rec(poly1); 1342 if (!rec) 1343 goto error; 1344 1345 for (i = 0; i < rec->n; ++i) { 1346 rec->p[i] = isl_poly_mul(rec->p[i], 1347 isl_poly_copy(poly2)); 1348 if (!rec->p[i]) 1349 goto error; 1350 } 1351 isl_poly_free(poly2); 1352 return poly1; 1353 } 1354 1355 is_cst = isl_poly_is_cst(poly1); 1356 if (is_cst < 0) 1357 goto error; 1358 if (is_cst) 1359 return isl_poly_mul_cst(poly1, poly2); 1360 1361 return isl_poly_mul_rec(poly1, poly2); 1362error: 1363 isl_poly_free(poly1); 1364 isl_poly_free(poly2); 1365 return NULL; 1366} 1367 1368__isl_give isl_poly *isl_poly_pow(__isl_take isl_poly *poly, unsigned power) 1369{ 1370 isl_poly *res; 1371 1372 if (!poly) 1373 return NULL; 1374 if (power == 1) 1375 return poly; 1376 1377 if (power % 2) 1378 res = isl_poly_copy(poly); 1379 else 1380 res = isl_poly_one(poly->ctx); 1381 1382 while (power >>= 1) { 1383 poly = isl_poly_mul(poly, isl_poly_copy(poly)); 1384 if (power % 2) 1385 res = isl_poly_mul(res, isl_poly_copy(poly)); 1386 } 1387 1388 isl_poly_free(poly); 1389 return res; 1390} 1391 1392__isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *space, 1393 unsigned n_div, __isl_take isl_poly *poly) 1394{ 1395 struct isl_qpolynomial *qp = NULL; 1396 isl_size total; 1397 1398 total = isl_space_dim(space, isl_dim_all); 1399 if (total < 0 || !poly) 1400 goto error; 1401 1402 if (!isl_space_is_set(space)) 1403 isl_die(isl_space_get_ctx(space), isl_error_invalid, 1404 "domain of polynomial should be a set", goto error); 1405 1406 qp = isl_calloc_type(space->ctx, struct isl_qpolynomial); 1407 if (!qp) 1408 goto error; 1409 1410 qp->ref = 1; 1411 qp->div = isl_mat_alloc(space->ctx, n_div, 1 + 1 + total + n_div); 1412 if (!qp->div) 1413 goto error; 1414 1415 qp->dim = space; 1416 qp->poly = poly; 1417 1418 return qp; 1419error: 1420 isl_space_free(space); 1421 isl_poly_free(poly); 1422 isl_qpolynomial_free(qp); 1423 return NULL; 1424} 1425 1426__isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp) 1427{ 1428 if (!qp) 1429 return NULL; 1430 1431 qp->ref++; 1432 return qp; 1433} 1434 1435/* Return a copy of the polynomial expression of "qp". 1436 */ 1437__isl_give isl_poly *isl_qpolynomial_get_poly(__isl_keep isl_qpolynomial *qp) 1438{ 1439 return qp ? isl_poly_copy(qp->poly) : NULL; 1440} 1441 1442/* Return the polynomial expression of "qp". 1443 * This may be either a copy or the polynomial expression itself 1444 * if there is only one reference to "qp". 1445 * This allows the polynomial expression to be modified inplace 1446 * if both the quasi-polynomial and its polynomial expression 1447 * have only a single reference. 1448 * The caller is not allowed to modify "qp" between this call and 1449 * a subsequent call to isl_qpolynomial_restore_poly. 1450 * The only exception is that isl_qpolynomial_free can be called instead. 1451 */ 1452static __isl_give isl_poly *isl_qpolynomial_take_poly( 1453 __isl_keep isl_qpolynomial *qp) 1454{ 1455 isl_poly *poly; 1456 1457 if (!qp) 1458 return NULL; 1459 if (qp->ref != 1) 1460 return isl_qpolynomial_get_poly(qp); 1461 poly = qp->poly; 1462 qp->poly = NULL; 1463 return poly; 1464} 1465 1466/* Set the polynomial expression of "qp" to "space", 1467 * where the polynomial expression of "qp" may be missing 1468 * due to a preceding call to isl_qpolynomial_take_poly. 1469 * However, in this case, "qp" only has a single reference and 1470 * then the call to isl_qpolynomial_cow has no effect. 1471 */ 1472static __isl_give isl_qpolynomial *isl_qpolynomial_restore_poly( 1473 __isl_keep isl_qpolynomial *qp, __isl_take isl_poly *poly) 1474{ 1475 if (!qp || !poly) 1476 goto error; 1477 1478 if (qp->poly == poly) { 1479 isl_poly_free(poly); 1480 return qp; 1481 } 1482 1483 qp = isl_qpolynomial_cow(qp); 1484 if (!qp) 1485 goto error; 1486 isl_poly_free(qp->poly); 1487 qp->poly = poly; 1488 1489 return qp; 1490error: 1491 isl_qpolynomial_free(qp); 1492 isl_poly_free(poly); 1493 return NULL; 1494} 1495 1496__isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp) 1497{ 1498 isl_poly *poly; 1499 struct isl_qpolynomial *dup; 1500 1501 if (!qp) 1502 return NULL; 1503 1504 poly = isl_qpolynomial_get_poly(qp); 1505 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, 1506 poly); 1507 if (!dup) 1508 return NULL; 1509 isl_mat_free(dup->div); 1510 dup->div = isl_qpolynomial_get_local(qp); 1511 if (!dup->div) 1512 goto error; 1513 1514 return dup; 1515error: 1516 isl_qpolynomial_free(dup); 1517 return NULL; 1518} 1519 1520__isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp) 1521{ 1522 if (!qp) 1523 return NULL; 1524 1525 if (qp->ref == 1) 1526 return qp; 1527 qp->ref--; 1528 return isl_qpolynomial_dup(qp); 1529} 1530 1531__isl_null isl_qpolynomial *isl_qpolynomial_free( 1532 __isl_take isl_qpolynomial *qp) 1533{ 1534 if (!qp) 1535 return NULL; 1536 1537 if (--qp->ref > 0) 1538 return NULL; 1539 1540 isl_space_free(qp->dim); 1541 isl_mat_free(qp->div); 1542 isl_poly_free(qp->poly); 1543 1544 free(qp); 1545 return NULL; 1546} 1547 1548__isl_give isl_poly *isl_poly_var_pow(isl_ctx *ctx, int pos, int power) 1549{ 1550 int i; 1551 isl_poly_rec *rec; 1552 isl_poly_cst *cst; 1553 1554 rec = isl_poly_alloc_rec(ctx, pos, 1 + power); 1555 if (!rec) 1556 return NULL; 1557 for (i = 0; i < 1 + power; ++i) { 1558 rec->p[i] = isl_poly_zero(ctx); 1559 if (!rec->p[i]) 1560 goto error; 1561 rec->n++; 1562 } 1563 cst = isl_poly_as_cst(rec->p[power]); 1564 isl_int_set_si(cst->n, 1); 1565 1566 return &rec->poly; 1567error: 1568 isl_poly_free(&rec->poly); 1569 return NULL; 1570} 1571 1572/* r array maps original positions to new positions. 1573 */ 1574static __isl_give isl_poly *reorder(__isl_take isl_poly *poly, int *r) 1575{ 1576 int i; 1577 isl_bool is_cst; 1578 isl_poly_rec *rec; 1579 isl_poly *base; 1580 isl_poly *res; 1581 1582 is_cst = isl_poly_is_cst(poly); 1583 if (is_cst < 0) 1584 return isl_poly_free(poly); 1585 if (is_cst) 1586 return poly; 1587 1588 rec = isl_poly_as_rec(poly); 1589 if (!rec) 1590 goto error; 1591 1592 isl_assert(poly->ctx, rec->n >= 1, goto error); 1593 1594 base = isl_poly_var_pow(poly->ctx, r[poly->var], 1); 1595 res = reorder(isl_poly_copy(rec->p[rec->n - 1]), r); 1596 1597 for (i = rec->n - 2; i >= 0; --i) { 1598 res = isl_poly_mul(res, isl_poly_copy(base)); 1599 res = isl_poly_sum(res, reorder(isl_poly_copy(rec->p[i]), r)); 1600 } 1601 1602 isl_poly_free(base); 1603 isl_poly_free(poly); 1604 1605 return res; 1606error: 1607 isl_poly_free(poly); 1608 return NULL; 1609} 1610 1611static isl_bool compatible_divs(__isl_keep isl_mat *div1, 1612 __isl_keep isl_mat *div2) 1613{ 1614 int n_row, n_col; 1615 isl_bool equal; 1616 1617 isl_assert(div1->ctx, div1->n_row >= div2->n_row && 1618 div1->n_col >= div2->n_col, 1619 return isl_bool_error); 1620 1621 if (div1->n_row == div2->n_row) 1622 return isl_mat_is_equal(div1, div2); 1623 1624 n_row = div1->n_row; 1625 n_col = div1->n_col; 1626 div1->n_row = div2->n_row; 1627 div1->n_col = div2->n_col; 1628 1629 equal = isl_mat_is_equal(div1, div2); 1630 1631 div1->n_row = n_row; 1632 div1->n_col = n_col; 1633 1634 return equal; 1635} 1636 1637static int cmp_row(__isl_keep isl_mat *div, int i, int j) 1638{ 1639 int li, lj; 1640 1641 li = isl_seq_last_non_zero(div->row[i], div->n_col); 1642 lj = isl_seq_last_non_zero(div->row[j], div->n_col); 1643 1644 if (li != lj) 1645 return li - lj; 1646 1647 return isl_seq_cmp(div->row[i], div->row[j], div->n_col); 1648} 1649 1650struct isl_div_sort_info { 1651 isl_mat *div; 1652 int row; 1653}; 1654 1655static int div_sort_cmp(const void *p1, const void *p2) 1656{ 1657 const struct isl_div_sort_info *i1, *i2; 1658 i1 = (const struct isl_div_sort_info *) p1; 1659 i2 = (const struct isl_div_sort_info *) p2; 1660 1661 return cmp_row(i1->div, i1->row, i2->row); 1662} 1663 1664/* Sort divs and remove duplicates. 1665 */ 1666static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp) 1667{ 1668 int i; 1669 int skip; 1670 int len; 1671 struct isl_div_sort_info *array = NULL; 1672 int *pos = NULL, *at = NULL; 1673 int *reordering = NULL; 1674 isl_size div_pos; 1675 1676 if (!qp) 1677 return NULL; 1678 if (qp->div->n_row <= 1) 1679 return qp; 1680 1681 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 1682 if (div_pos < 0) 1683 return isl_qpolynomial_free(qp); 1684 1685 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info, 1686 qp->div->n_row); 1687 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row); 1688 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row); 1689 len = qp->div->n_col - 2; 1690 reordering = isl_alloc_array(qp->div->ctx, int, len); 1691 if (!array || !pos || !at || !reordering) 1692 goto error; 1693 1694 for (i = 0; i < qp->div->n_row; ++i) { 1695 array[i].div = qp->div; 1696 array[i].row = i; 1697 pos[i] = i; 1698 at[i] = i; 1699 } 1700 1701 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info), 1702 div_sort_cmp); 1703 1704 for (i = 0; i < div_pos; ++i) 1705 reordering[i] = i; 1706 1707 for (i = 0; i < qp->div->n_row; ++i) { 1708 if (pos[array[i].row] == i) 1709 continue; 1710 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]); 1711 pos[at[i]] = pos[array[i].row]; 1712 at[pos[array[i].row]] = at[i]; 1713 at[i] = array[i].row; 1714 pos[array[i].row] = i; 1715 } 1716 1717 skip = 0; 1718 for (i = 0; i < len - div_pos; ++i) { 1719 if (i > 0 && 1720 isl_seq_eq(qp->div->row[i - skip - 1], 1721 qp->div->row[i - skip], qp->div->n_col)) { 1722 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1); 1723 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1, 1724 2 + div_pos + i - skip); 1725 qp->div = isl_mat_drop_cols(qp->div, 1726 2 + div_pos + i - skip, 1); 1727 skip++; 1728 } 1729 reordering[div_pos + array[i].row] = div_pos + i - skip; 1730 } 1731 1732 qp->poly = reorder(qp->poly, reordering); 1733 1734 if (!qp->poly || !qp->div) 1735 goto error; 1736 1737 free(at); 1738 free(pos); 1739 free(array); 1740 free(reordering); 1741 1742 return qp; 1743error: 1744 free(at); 1745 free(pos); 1746 free(array); 1747 free(reordering); 1748 isl_qpolynomial_free(qp); 1749 return NULL; 1750} 1751 1752static __isl_give isl_poly *expand(__isl_take isl_poly *poly, int *exp, 1753 int first) 1754{ 1755 int i; 1756 isl_bool is_cst; 1757 isl_poly_rec *rec; 1758 1759 is_cst = isl_poly_is_cst(poly); 1760 if (is_cst < 0) 1761 return isl_poly_free(poly); 1762 if (is_cst) 1763 return poly; 1764 1765 if (poly->var < first) 1766 return poly; 1767 1768 if (exp[poly->var - first] == poly->var - first) 1769 return poly; 1770 1771 poly = isl_poly_cow(poly); 1772 if (!poly) 1773 goto error; 1774 1775 poly->var = exp[poly->var - first] + first; 1776 1777 rec = isl_poly_as_rec(poly); 1778 if (!rec) 1779 goto error; 1780 1781 for (i = 0; i < rec->n; ++i) { 1782 rec->p[i] = expand(rec->p[i], exp, first); 1783 if (!rec->p[i]) 1784 goto error; 1785 } 1786 1787 return poly; 1788error: 1789 isl_poly_free(poly); 1790 return NULL; 1791} 1792 1793static __isl_give isl_qpolynomial *with_merged_divs( 1794 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1, 1795 __isl_take isl_qpolynomial *qp2), 1796 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2) 1797{ 1798 int *exp1 = NULL; 1799 int *exp2 = NULL; 1800 isl_mat *div = NULL; 1801 int n_div1, n_div2; 1802 1803 qp1 = isl_qpolynomial_cow(qp1); 1804 qp2 = isl_qpolynomial_cow(qp2); 1805 1806 if (!qp1 || !qp2) 1807 goto error; 1808 1809 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row && 1810 qp1->div->n_col >= qp2->div->n_col, goto error); 1811 1812 n_div1 = qp1->div->n_row; 1813 n_div2 = qp2->div->n_row; 1814 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1); 1815 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2); 1816 if ((n_div1 && !exp1) || (n_div2 && !exp2)) 1817 goto error; 1818 1819 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2); 1820 if (!div) 1821 goto error; 1822 1823 isl_mat_free(qp1->div); 1824 qp1->div = isl_mat_copy(div); 1825 isl_mat_free(qp2->div); 1826 qp2->div = isl_mat_copy(div); 1827 1828 qp1->poly = expand(qp1->poly, exp1, div->n_col - div->n_row - 2); 1829 qp2->poly = expand(qp2->poly, exp2, div->n_col - div->n_row - 2); 1830 1831 if (!qp1->poly || !qp2->poly) 1832 goto error; 1833 1834 isl_mat_free(div); 1835 free(exp1); 1836 free(exp2); 1837 1838 return fn(qp1, qp2); 1839error: 1840 isl_mat_free(div); 1841 free(exp1); 1842 free(exp2); 1843 isl_qpolynomial_free(qp1); 1844 isl_qpolynomial_free(qp2); 1845 return NULL; 1846} 1847 1848__isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1, 1849 __isl_take isl_qpolynomial *qp2) 1850{ 1851 isl_bool compatible; 1852 isl_poly *poly; 1853 1854 if (isl_qpolynomial_check_equal_space(qp1, qp2) < 0) 1855 goto error; 1856 1857 if (qp1->div->n_row < qp2->div->n_row) 1858 return isl_qpolynomial_add(qp2, qp1); 1859 1860 compatible = compatible_divs(qp1->div, qp2->div); 1861 if (compatible < 0) 1862 goto error; 1863 if (!compatible) 1864 return with_merged_divs(isl_qpolynomial_add, qp1, qp2); 1865 1866 poly = isl_qpolynomial_take_poly(qp1); 1867 poly = isl_poly_sum(poly, isl_qpolynomial_get_poly(qp2)); 1868 qp1 = isl_qpolynomial_restore_poly(qp1, poly); 1869 1870 isl_qpolynomial_free(qp2); 1871 1872 return qp1; 1873error: 1874 isl_qpolynomial_free(qp1); 1875 isl_qpolynomial_free(qp2); 1876 return NULL; 1877} 1878 1879__isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain( 1880 __isl_keep isl_set *dom, 1881 __isl_take isl_qpolynomial *qp1, 1882 __isl_take isl_qpolynomial *qp2) 1883{ 1884 qp1 = isl_qpolynomial_add(qp1, qp2); 1885 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom)); 1886 return qp1; 1887} 1888 1889__isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1, 1890 __isl_take isl_qpolynomial *qp2) 1891{ 1892 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2)); 1893} 1894 1895__isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int( 1896 __isl_take isl_qpolynomial *qp, isl_int v) 1897{ 1898 isl_poly *poly; 1899 1900 if (isl_int_is_zero(v)) 1901 return qp; 1902 1903 poly = isl_qpolynomial_take_poly(qp); 1904 poly = isl_poly_add_isl_int(poly, v); 1905 qp = isl_qpolynomial_restore_poly(qp, poly); 1906 1907 return qp; 1908} 1909 1910__isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp) 1911{ 1912 if (!qp) 1913 return NULL; 1914 1915 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone); 1916} 1917 1918__isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int( 1919 __isl_take isl_qpolynomial *qp, isl_int v) 1920{ 1921 isl_poly *poly; 1922 1923 if (isl_int_is_one(v)) 1924 return qp; 1925 1926 if (qp && isl_int_is_zero(v)) { 1927 isl_qpolynomial *zero; 1928 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim)); 1929 isl_qpolynomial_free(qp); 1930 return zero; 1931 } 1932 1933 poly = isl_qpolynomial_take_poly(qp); 1934 poly = isl_poly_mul_isl_int(poly, v); 1935 qp = isl_qpolynomial_restore_poly(qp, poly); 1936 1937 return qp; 1938} 1939 1940__isl_give isl_qpolynomial *isl_qpolynomial_scale( 1941 __isl_take isl_qpolynomial *qp, isl_int v) 1942{ 1943 return isl_qpolynomial_mul_isl_int(qp, v); 1944} 1945 1946/* Multiply "qp" by "v". 1947 */ 1948__isl_give isl_qpolynomial *isl_qpolynomial_scale_val( 1949 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v) 1950{ 1951 isl_poly *poly; 1952 1953 if (!qp || !v) 1954 goto error; 1955 1956 if (!isl_val_is_rat(v)) 1957 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid, 1958 "expecting rational factor", goto error); 1959 1960 if (isl_val_is_one(v)) { 1961 isl_val_free(v); 1962 return qp; 1963 } 1964 1965 if (isl_val_is_zero(v)) { 1966 isl_space *space; 1967 1968 space = isl_qpolynomial_get_domain_space(qp); 1969 isl_qpolynomial_free(qp); 1970 isl_val_free(v); 1971 return isl_qpolynomial_zero_on_domain(space); 1972 } 1973 1974 poly = isl_qpolynomial_take_poly(qp); 1975 poly = isl_poly_scale_val(poly, v); 1976 qp = isl_qpolynomial_restore_poly(qp, poly); 1977 1978 isl_val_free(v); 1979 return qp; 1980error: 1981 isl_val_free(v); 1982 isl_qpolynomial_free(qp); 1983 return NULL; 1984} 1985 1986/* Divide "qp" by "v". 1987 */ 1988__isl_give isl_qpolynomial *isl_qpolynomial_scale_down_val( 1989 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v) 1990{ 1991 if (!qp || !v) 1992 goto error; 1993 1994 if (!isl_val_is_rat(v)) 1995 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid, 1996 "expecting rational factor", goto error); 1997 if (isl_val_is_zero(v)) 1998 isl_die(isl_val_get_ctx(v), isl_error_invalid, 1999 "cannot scale down by zero", goto error); 2000 2001 return isl_qpolynomial_scale_val(qp, isl_val_inv(v)); 2002error: 2003 isl_val_free(v); 2004 isl_qpolynomial_free(qp); 2005 return NULL; 2006} 2007 2008__isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1, 2009 __isl_take isl_qpolynomial *qp2) 2010{ 2011 isl_bool compatible; 2012 isl_poly *poly; 2013 2014 if (isl_qpolynomial_check_equal_space(qp1, qp2) < 0) 2015 goto error; 2016 2017 if (qp1->div->n_row < qp2->div->n_row) 2018 return isl_qpolynomial_mul(qp2, qp1); 2019 2020 compatible = compatible_divs(qp1->div, qp2->div); 2021 if (compatible < 0) 2022 goto error; 2023 if (!compatible) 2024 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2); 2025 2026 poly = isl_qpolynomial_take_poly(qp1); 2027 poly = isl_poly_mul(poly, isl_qpolynomial_get_poly(qp2)); 2028 qp1 = isl_qpolynomial_restore_poly(qp1, poly); 2029 2030 isl_qpolynomial_free(qp2); 2031 2032 return qp1; 2033error: 2034 isl_qpolynomial_free(qp1); 2035 isl_qpolynomial_free(qp2); 2036 return NULL; 2037} 2038 2039__isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp, 2040 unsigned power) 2041{ 2042 isl_poly *poly; 2043 2044 poly = isl_qpolynomial_take_poly(qp); 2045 poly = isl_poly_pow(poly, power); 2046 qp = isl_qpolynomial_restore_poly(qp, poly); 2047 2048 return qp; 2049} 2050 2051__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow( 2052 __isl_take isl_pw_qpolynomial *pwqp, unsigned power) 2053{ 2054 int i; 2055 2056 if (power == 1) 2057 return pwqp; 2058 2059 pwqp = isl_pw_qpolynomial_cow(pwqp); 2060 if (!pwqp) 2061 return NULL; 2062 2063 for (i = 0; i < pwqp->n; ++i) { 2064 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power); 2065 if (!pwqp->p[i].qp) 2066 return isl_pw_qpolynomial_free(pwqp); 2067 } 2068 2069 return pwqp; 2070} 2071 2072__isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain( 2073 __isl_take isl_space *domain) 2074{ 2075 if (!domain) 2076 return NULL; 2077 return isl_qpolynomial_alloc(domain, 0, isl_poly_zero(domain->ctx)); 2078} 2079 2080__isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain( 2081 __isl_take isl_space *domain) 2082{ 2083 if (!domain) 2084 return NULL; 2085 return isl_qpolynomial_alloc(domain, 0, isl_poly_one(domain->ctx)); 2086} 2087 2088__isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain( 2089 __isl_take isl_space *domain) 2090{ 2091 if (!domain) 2092 return NULL; 2093 return isl_qpolynomial_alloc(domain, 0, isl_poly_infty(domain->ctx)); 2094} 2095 2096__isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain( 2097 __isl_take isl_space *domain) 2098{ 2099 if (!domain) 2100 return NULL; 2101 return isl_qpolynomial_alloc(domain, 0, isl_poly_neginfty(domain->ctx)); 2102} 2103 2104__isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain( 2105 __isl_take isl_space *domain) 2106{ 2107 if (!domain) 2108 return NULL; 2109 return isl_qpolynomial_alloc(domain, 0, isl_poly_nan(domain->ctx)); 2110} 2111 2112__isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain( 2113 __isl_take isl_space *domain, 2114 isl_int v) 2115{ 2116 struct isl_qpolynomial *qp; 2117 isl_poly_cst *cst; 2118 2119 qp = isl_qpolynomial_zero_on_domain(domain); 2120 if (!qp) 2121 return NULL; 2122 2123 cst = isl_poly_as_cst(qp->poly); 2124 isl_int_set(cst->n, v); 2125 2126 return qp; 2127} 2128 2129isl_bool isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp, 2130 isl_int *n, isl_int *d) 2131{ 2132 isl_bool is_cst; 2133 isl_poly_cst *cst; 2134 2135 if (!qp) 2136 return isl_bool_error; 2137 2138 is_cst = isl_poly_is_cst(qp->poly); 2139 if (is_cst < 0 || !is_cst) 2140 return is_cst; 2141 2142 cst = isl_poly_as_cst(qp->poly); 2143 if (!cst) 2144 return isl_bool_error; 2145 2146 if (n) 2147 isl_int_set(*n, cst->n); 2148 if (d) 2149 isl_int_set(*d, cst->d); 2150 2151 return isl_bool_true; 2152} 2153 2154/* Return the constant term of "poly". 2155 */ 2156static __isl_give isl_val *isl_poly_get_constant_val(__isl_keep isl_poly *poly) 2157{ 2158 isl_bool is_cst; 2159 isl_poly_cst *cst; 2160 2161 if (!poly) 2162 return NULL; 2163 2164 while ((is_cst = isl_poly_is_cst(poly)) == isl_bool_false) { 2165 isl_poly_rec *rec; 2166 2167 rec = isl_poly_as_rec(poly); 2168 if (!rec) 2169 return NULL; 2170 poly = rec->p[0]; 2171 } 2172 if (is_cst < 0) 2173 return NULL; 2174 2175 cst = isl_poly_as_cst(poly); 2176 if (!cst) 2177 return NULL; 2178 return isl_val_rat_from_isl_int(cst->poly.ctx, cst->n, cst->d); 2179} 2180 2181/* Return the constant term of "qp". 2182 */ 2183__isl_give isl_val *isl_qpolynomial_get_constant_val( 2184 __isl_keep isl_qpolynomial *qp) 2185{ 2186 if (!qp) 2187 return NULL; 2188 2189 return isl_poly_get_constant_val(qp->poly); 2190} 2191 2192isl_bool isl_poly_is_affine(__isl_keep isl_poly *poly) 2193{ 2194 isl_bool is_cst; 2195 isl_poly_rec *rec; 2196 2197 if (!poly) 2198 return isl_bool_error; 2199 2200 if (poly->var < 0) 2201 return isl_bool_true; 2202 2203 rec = isl_poly_as_rec(poly); 2204 if (!rec) 2205 return isl_bool_error; 2206 2207 if (rec->n > 2) 2208 return isl_bool_false; 2209 2210 isl_assert(poly->ctx, rec->n > 1, return isl_bool_error); 2211 2212 is_cst = isl_poly_is_cst(rec->p[1]); 2213 if (is_cst < 0 || !is_cst) 2214 return is_cst; 2215 2216 return isl_poly_is_affine(rec->p[0]); 2217} 2218 2219isl_bool isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp) 2220{ 2221 if (!qp) 2222 return isl_bool_error; 2223 2224 if (qp->div->n_row > 0) 2225 return isl_bool_false; 2226 2227 return isl_poly_is_affine(qp->poly); 2228} 2229 2230static void update_coeff(__isl_keep isl_vec *aff, 2231 __isl_keep isl_poly_cst *cst, int pos) 2232{ 2233 isl_int gcd; 2234 isl_int f; 2235 2236 if (isl_int_is_zero(cst->n)) 2237 return; 2238 2239 isl_int_init(gcd); 2240 isl_int_init(f); 2241 isl_int_gcd(gcd, cst->d, aff->el[0]); 2242 isl_int_divexact(f, cst->d, gcd); 2243 isl_int_divexact(gcd, aff->el[0], gcd); 2244 isl_seq_scale(aff->el, aff->el, f, aff->size); 2245 isl_int_mul(aff->el[1 + pos], gcd, cst->n); 2246 isl_int_clear(gcd); 2247 isl_int_clear(f); 2248} 2249 2250int isl_poly_update_affine(__isl_keep isl_poly *poly, __isl_keep isl_vec *aff) 2251{ 2252 isl_poly_cst *cst; 2253 isl_poly_rec *rec; 2254 2255 if (!poly || !aff) 2256 return -1; 2257 2258 if (poly->var < 0) { 2259 isl_poly_cst *cst; 2260 2261 cst = isl_poly_as_cst(poly); 2262 if (!cst) 2263 return -1; 2264 update_coeff(aff, cst, 0); 2265 return 0; 2266 } 2267 2268 rec = isl_poly_as_rec(poly); 2269 if (!rec) 2270 return -1; 2271 isl_assert(poly->ctx, rec->n == 2, return -1); 2272 2273 cst = isl_poly_as_cst(rec->p[1]); 2274 if (!cst) 2275 return -1; 2276 update_coeff(aff, cst, 1 + poly->var); 2277 2278 return isl_poly_update_affine(rec->p[0], aff); 2279} 2280 2281__isl_give isl_vec *isl_qpolynomial_extract_affine( 2282 __isl_keep isl_qpolynomial *qp) 2283{ 2284 isl_vec *aff; 2285 isl_size d; 2286 2287 d = isl_qpolynomial_domain_dim(qp, isl_dim_all); 2288 if (d < 0) 2289 return NULL; 2290 2291 aff = isl_vec_alloc(qp->div->ctx, 2 + d); 2292 if (!aff) 2293 return NULL; 2294 2295 isl_seq_clr(aff->el + 1, 1 + d); 2296 isl_int_set_si(aff->el[0], 1); 2297 2298 if (isl_poly_update_affine(qp->poly, aff) < 0) 2299 goto error; 2300 2301 return aff; 2302error: 2303 isl_vec_free(aff); 2304 return NULL; 2305} 2306 2307/* Compare two quasi-polynomials. 2308 * 2309 * Return -1 if "qp1" is "smaller" than "qp2", 1 if "qp1" is "greater" 2310 * than "qp2" and 0 if they are equal. 2311 */ 2312int isl_qpolynomial_plain_cmp(__isl_keep isl_qpolynomial *qp1, 2313 __isl_keep isl_qpolynomial *qp2) 2314{ 2315 int cmp; 2316 2317 if (qp1 == qp2) 2318 return 0; 2319 if (!qp1) 2320 return -1; 2321 if (!qp2) 2322 return 1; 2323 2324 cmp = isl_space_cmp(qp1->dim, qp2->dim); 2325 if (cmp != 0) 2326 return cmp; 2327 2328 cmp = isl_local_cmp(qp1->div, qp2->div); 2329 if (cmp != 0) 2330 return cmp; 2331 2332 return isl_poly_plain_cmp(qp1->poly, qp2->poly); 2333} 2334 2335/* Is "qp1" obviously equal to "qp2"? 2336 * 2337 * NaN is not equal to anything, not even to another NaN. 2338 */ 2339isl_bool isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1, 2340 __isl_keep isl_qpolynomial *qp2) 2341{ 2342 isl_bool equal; 2343 2344 if (!qp1 || !qp2) 2345 return isl_bool_error; 2346 2347 if (isl_qpolynomial_is_nan(qp1) || isl_qpolynomial_is_nan(qp2)) 2348 return isl_bool_false; 2349 2350 equal = isl_space_is_equal(qp1->dim, qp2->dim); 2351 if (equal < 0 || !equal) 2352 return equal; 2353 2354 equal = isl_mat_is_equal(qp1->div, qp2->div); 2355 if (equal < 0 || !equal) 2356 return equal; 2357 2358 return isl_poly_is_equal(qp1->poly, qp2->poly); 2359} 2360 2361static isl_stat poly_update_den(__isl_keep isl_poly *poly, isl_int *d) 2362{ 2363 int i; 2364 isl_bool is_cst; 2365 isl_poly_rec *rec; 2366 2367 is_cst = isl_poly_is_cst(poly); 2368 if (is_cst < 0) 2369 return isl_stat_error; 2370 if (is_cst) { 2371 isl_poly_cst *cst; 2372 cst = isl_poly_as_cst(poly); 2373 if (!cst) 2374 return isl_stat_error; 2375 isl_int_lcm(*d, *d, cst->d); 2376 return isl_stat_ok; 2377 } 2378 2379 rec = isl_poly_as_rec(poly); 2380 if (!rec) 2381 return isl_stat_error; 2382 2383 for (i = 0; i < rec->n; ++i) 2384 poly_update_den(rec->p[i], d); 2385 2386 return isl_stat_ok; 2387} 2388 2389__isl_give isl_val *isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp) 2390{ 2391 isl_val *d; 2392 2393 if (!qp) 2394 return NULL; 2395 d = isl_val_one(isl_qpolynomial_get_ctx(qp)); 2396 if (!d) 2397 return NULL; 2398 if (poly_update_den(qp->poly, &d->n) < 0) 2399 return isl_val_free(d); 2400 return d; 2401} 2402 2403__isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain( 2404 __isl_take isl_space *domain, int pos, int power) 2405{ 2406 struct isl_ctx *ctx; 2407 2408 if (!domain) 2409 return NULL; 2410 2411 ctx = domain->ctx; 2412 2413 return isl_qpolynomial_alloc(domain, 0, 2414 isl_poly_var_pow(ctx, pos, power)); 2415} 2416 2417__isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain( 2418 __isl_take isl_space *domain, enum isl_dim_type type, unsigned pos) 2419{ 2420 isl_size off; 2421 2422 if (isl_space_check_is_set(domain ) < 0) 2423 goto error; 2424 if (isl_space_check_range(domain, type, pos, 1) < 0) 2425 goto error; 2426 2427 off = isl_space_offset(domain, type); 2428 if (off < 0) 2429 goto error; 2430 2431 return isl_qpolynomial_var_pow_on_domain(domain, off + pos, 1); 2432error: 2433 isl_space_free(domain); 2434 return NULL; 2435} 2436 2437__isl_give isl_poly *isl_poly_subs(__isl_take isl_poly *poly, 2438 unsigned first, unsigned n, __isl_keep isl_poly **subs) 2439{ 2440 int i; 2441 isl_bool is_cst; 2442 isl_poly_rec *rec; 2443 isl_poly *base, *res; 2444 2445 is_cst = isl_poly_is_cst(poly); 2446 if (is_cst < 0) 2447 return isl_poly_free(poly); 2448 if (is_cst) 2449 return poly; 2450 2451 if (poly->var < first) 2452 return poly; 2453 2454 rec = isl_poly_as_rec(poly); 2455 if (!rec) 2456 goto error; 2457 2458 isl_assert(poly->ctx, rec->n >= 1, goto error); 2459 2460 if (poly->var >= first + n) 2461 base = isl_poly_var_pow(poly->ctx, poly->var, 1); 2462 else 2463 base = isl_poly_copy(subs[poly->var - first]); 2464 2465 res = isl_poly_subs(isl_poly_copy(rec->p[rec->n - 1]), first, n, subs); 2466 for (i = rec->n - 2; i >= 0; --i) { 2467 isl_poly *t; 2468 t = isl_poly_subs(isl_poly_copy(rec->p[i]), first, n, subs); 2469 res = isl_poly_mul(res, isl_poly_copy(base)); 2470 res = isl_poly_sum(res, t); 2471 } 2472 2473 isl_poly_free(base); 2474 isl_poly_free(poly); 2475 2476 return res; 2477error: 2478 isl_poly_free(poly); 2479 return NULL; 2480} 2481 2482__isl_give isl_poly *isl_poly_from_affine(isl_ctx *ctx, isl_int *f, 2483 isl_int denom, unsigned len) 2484{ 2485 int i; 2486 isl_poly *poly; 2487 2488 isl_assert(ctx, len >= 1, return NULL); 2489 2490 poly = isl_poly_rat_cst(ctx, f[0], denom); 2491 for (i = 0; i < len - 1; ++i) { 2492 isl_poly *t; 2493 isl_poly *c; 2494 2495 if (isl_int_is_zero(f[1 + i])) 2496 continue; 2497 2498 c = isl_poly_rat_cst(ctx, f[1 + i], denom); 2499 t = isl_poly_var_pow(ctx, i, 1); 2500 t = isl_poly_mul(c, t); 2501 poly = isl_poly_sum(poly, t); 2502 } 2503 2504 return poly; 2505} 2506 2507/* Remove common factor of non-constant terms and denominator. 2508 */ 2509static void normalize_div(__isl_keep isl_qpolynomial *qp, int div) 2510{ 2511 isl_ctx *ctx = qp->div->ctx; 2512 unsigned total = qp->div->n_col - 2; 2513 2514 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd); 2515 isl_int_gcd(ctx->normalize_gcd, 2516 ctx->normalize_gcd, qp->div->row[div][0]); 2517 if (isl_int_is_one(ctx->normalize_gcd)) 2518 return; 2519 2520 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2, 2521 ctx->normalize_gcd, total); 2522 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0], 2523 ctx->normalize_gcd); 2524 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1], 2525 ctx->normalize_gcd); 2526} 2527 2528/* Replace the integer division identified by "div" by the polynomial "s". 2529 * The integer division is assumed not to appear in the definition 2530 * of any other integer divisions. 2531 */ 2532static __isl_give isl_qpolynomial *substitute_div( 2533 __isl_take isl_qpolynomial *qp, int div, __isl_take isl_poly *s) 2534{ 2535 int i; 2536 isl_size div_pos; 2537 int *reordering; 2538 isl_ctx *ctx; 2539 2540 if (!qp || !s) 2541 goto error; 2542 2543 qp = isl_qpolynomial_cow(qp); 2544 if (!qp) 2545 goto error; 2546 2547 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 2548 if (div_pos < 0) 2549 goto error; 2550 qp->poly = isl_poly_subs(qp->poly, div_pos + div, 1, &s); 2551 if (!qp->poly) 2552 goto error; 2553 2554 ctx = isl_qpolynomial_get_ctx(qp); 2555 reordering = isl_alloc_array(ctx, int, div_pos + qp->div->n_row); 2556 if (!reordering) 2557 goto error; 2558 for (i = 0; i < div_pos + div; ++i) 2559 reordering[i] = i; 2560 for (i = div_pos + div + 1; i < div_pos + qp->div->n_row; ++i) 2561 reordering[i] = i - 1; 2562 qp->div = isl_mat_drop_rows(qp->div, div, 1); 2563 qp->div = isl_mat_drop_cols(qp->div, 2 + div_pos + div, 1); 2564 qp->poly = reorder(qp->poly, reordering); 2565 free(reordering); 2566 2567 if (!qp->poly || !qp->div) 2568 goto error; 2569 2570 isl_poly_free(s); 2571 return qp; 2572error: 2573 isl_qpolynomial_free(qp); 2574 isl_poly_free(s); 2575 return NULL; 2576} 2577 2578/* Replace all integer divisions [e/d] that turn out to not actually be integer 2579 * divisions because d is equal to 1 by their definition, i.e., e. 2580 */ 2581static __isl_give isl_qpolynomial *substitute_non_divs( 2582 __isl_take isl_qpolynomial *qp) 2583{ 2584 int i, j; 2585 isl_size div_pos; 2586 isl_poly *s; 2587 2588 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 2589 if (div_pos < 0) 2590 return isl_qpolynomial_free(qp); 2591 2592 for (i = 0; qp && i < qp->div->n_row; ++i) { 2593 if (!isl_int_is_one(qp->div->row[i][0])) 2594 continue; 2595 for (j = i + 1; j < qp->div->n_row; ++j) { 2596 if (isl_int_is_zero(qp->div->row[j][2 + div_pos + i])) 2597 continue; 2598 isl_seq_combine(qp->div->row[j] + 1, 2599 qp->div->ctx->one, qp->div->row[j] + 1, 2600 qp->div->row[j][2 + div_pos + i], 2601 qp->div->row[i] + 1, 1 + div_pos + i); 2602 isl_int_set_si(qp->div->row[j][2 + div_pos + i], 0); 2603 normalize_div(qp, j); 2604 } 2605 s = isl_poly_from_affine(qp->dim->ctx, qp->div->row[i] + 1, 2606 qp->div->row[i][0], qp->div->n_col - 1); 2607 qp = substitute_div(qp, i, s); 2608 --i; 2609 } 2610 2611 return qp; 2612} 2613 2614/* Reduce the coefficients of div "div" to lie in the interval [0, d-1], 2615 * with d the denominator. When replacing the coefficient e of x by 2616 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x 2617 * inside the division, so we need to add floor(e/d) * x outside. 2618 * That is, we replace q by q' + floor(e/d) * x and we therefore need 2619 * to adjust the coefficient of x in each later div that depends on the 2620 * current div "div" and also in the affine expressions in the rows of "mat" 2621 * (if they too depend on "div"). 2622 */ 2623static void reduce_div(__isl_keep isl_qpolynomial *qp, int div, 2624 __isl_keep isl_mat **mat) 2625{ 2626 int i, j; 2627 isl_int v; 2628 unsigned total = qp->div->n_col - qp->div->n_row - 2; 2629 2630 isl_int_init(v); 2631 for (i = 0; i < 1 + total + div; ++i) { 2632 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) && 2633 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0])) 2634 continue; 2635 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]); 2636 isl_int_fdiv_r(qp->div->row[div][1 + i], 2637 qp->div->row[div][1 + i], qp->div->row[div][0]); 2638 *mat = isl_mat_col_addmul(*mat, i, v, 1 + total + div); 2639 for (j = div + 1; j < qp->div->n_row; ++j) { 2640 if (isl_int_is_zero(qp->div->row[j][2 + total + div])) 2641 continue; 2642 isl_int_addmul(qp->div->row[j][1 + i], 2643 v, qp->div->row[j][2 + total + div]); 2644 } 2645 } 2646 isl_int_clear(v); 2647} 2648 2649/* Check if the last non-zero coefficient is bigger that half of the 2650 * denominator. If so, we will invert the div to further reduce the number 2651 * of distinct divs that may appear. 2652 * If the last non-zero coefficient is exactly half the denominator, 2653 * then we continue looking for earlier coefficients that are bigger 2654 * than half the denominator. 2655 */ 2656static int needs_invert(__isl_keep isl_mat *div, int row) 2657{ 2658 int i; 2659 int cmp; 2660 2661 for (i = div->n_col - 1; i >= 1; --i) { 2662 if (isl_int_is_zero(div->row[row][i])) 2663 continue; 2664 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2); 2665 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]); 2666 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2); 2667 if (cmp) 2668 return cmp > 0; 2669 if (i == 1) 2670 return 1; 2671 } 2672 2673 return 0; 2674} 2675 2676/* Replace div "div" q = [e/d] by -[(-e+(d-1))/d]. 2677 * We only invert the coefficients of e (and the coefficient of q in 2678 * later divs and in the rows of "mat"). After calling this function, the 2679 * coefficients of e should be reduced again. 2680 */ 2681static void invert_div(__isl_keep isl_qpolynomial *qp, int div, 2682 __isl_keep isl_mat **mat) 2683{ 2684 unsigned total = qp->div->n_col - qp->div->n_row - 2; 2685 2686 isl_seq_neg(qp->div->row[div] + 1, 2687 qp->div->row[div] + 1, qp->div->n_col - 1); 2688 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1); 2689 isl_int_add(qp->div->row[div][1], 2690 qp->div->row[div][1], qp->div->row[div][0]); 2691 *mat = isl_mat_col_neg(*mat, 1 + total + div); 2692 isl_mat_col_mul(qp->div, 2 + total + div, 2693 qp->div->ctx->negone, 2 + total + div); 2694} 2695 2696/* Reduce all divs of "qp" to have coefficients 2697 * in the interval [0, d-1], with d the denominator and such that the 2698 * last non-zero coefficient that is not equal to d/2 is smaller than d/2. 2699 * The modifications to the integer divisions need to be reflected 2700 * in the factors of the polynomial that refer to the original 2701 * integer divisions. To this end, the modifications are collected 2702 * as a set of affine expressions and then plugged into the polynomial. 2703 * 2704 * After the reduction, some divs may have become redundant or identical, 2705 * so we call substitute_non_divs and sort_divs. If these functions 2706 * eliminate divs or merge two or more divs into one, the coefficients 2707 * of the enclosing divs may have to be reduced again, so we call 2708 * ourselves recursively if the number of divs decreases. 2709 */ 2710static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp) 2711{ 2712 int i; 2713 isl_ctx *ctx; 2714 isl_mat *mat; 2715 isl_poly **s; 2716 unsigned o_div; 2717 isl_size n_div, total, new_n_div; 2718 2719 total = isl_qpolynomial_domain_dim(qp, isl_dim_all); 2720 n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div); 2721 o_div = isl_qpolynomial_domain_offset(qp, isl_dim_div); 2722 if (total < 0 || n_div < 0) 2723 return isl_qpolynomial_free(qp); 2724 ctx = isl_qpolynomial_get_ctx(qp); 2725 mat = isl_mat_zero(ctx, n_div, 1 + total); 2726 2727 for (i = 0; i < n_div; ++i) 2728 mat = isl_mat_set_element_si(mat, i, o_div + i, 1); 2729 2730 for (i = 0; i < qp->div->n_row; ++i) { 2731 normalize_div(qp, i); 2732 reduce_div(qp, i, &mat); 2733 if (needs_invert(qp->div, i)) { 2734 invert_div(qp, i, &mat); 2735 reduce_div(qp, i, &mat); 2736 } 2737 } 2738 if (!mat) 2739 goto error; 2740 2741 s = isl_alloc_array(ctx, struct isl_poly *, n_div); 2742 if (n_div && !s) 2743 goto error; 2744 for (i = 0; i < n_div; ++i) 2745 s[i] = isl_poly_from_affine(ctx, mat->row[i], ctx->one, 2746 1 + total); 2747 qp->poly = isl_poly_subs(qp->poly, o_div - 1, n_div, s); 2748 for (i = 0; i < n_div; ++i) 2749 isl_poly_free(s[i]); 2750 free(s); 2751 if (!qp->poly) 2752 goto error; 2753 2754 isl_mat_free(mat); 2755 2756 qp = substitute_non_divs(qp); 2757 qp = sort_divs(qp); 2758 new_n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div); 2759 if (new_n_div < 0) 2760 return isl_qpolynomial_free(qp); 2761 if (new_n_div < n_div) 2762 return reduce_divs(qp); 2763 2764 return qp; 2765error: 2766 isl_qpolynomial_free(qp); 2767 isl_mat_free(mat); 2768 return NULL; 2769} 2770 2771__isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain( 2772 __isl_take isl_space *domain, const isl_int n, const isl_int d) 2773{ 2774 struct isl_qpolynomial *qp; 2775 isl_poly_cst *cst; 2776 2777 qp = isl_qpolynomial_zero_on_domain(domain); 2778 if (!qp) 2779 return NULL; 2780 2781 cst = isl_poly_as_cst(qp->poly); 2782 isl_int_set(cst->n, n); 2783 isl_int_set(cst->d, d); 2784 2785 return qp; 2786} 2787 2788/* Return an isl_qpolynomial that is equal to "val" on domain space "domain". 2789 */ 2790__isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain( 2791 __isl_take isl_space *domain, __isl_take isl_val *val) 2792{ 2793 isl_qpolynomial *qp; 2794 isl_poly_cst *cst; 2795 2796 qp = isl_qpolynomial_zero_on_domain(domain); 2797 if (!qp || !val) 2798 goto error; 2799 2800 cst = isl_poly_as_cst(qp->poly); 2801 isl_int_set(cst->n, val->n); 2802 isl_int_set(cst->d, val->d); 2803 2804 isl_val_free(val); 2805 return qp; 2806error: 2807 isl_val_free(val); 2808 isl_qpolynomial_free(qp); 2809 return NULL; 2810} 2811 2812static isl_stat poly_set_active(__isl_keep isl_poly *poly, int *active, int d) 2813{ 2814 isl_bool is_cst; 2815 isl_poly_rec *rec; 2816 int i; 2817 2818 is_cst = isl_poly_is_cst(poly); 2819 if (is_cst < 0) 2820 return isl_stat_error; 2821 if (is_cst) 2822 return isl_stat_ok; 2823 2824 if (poly->var < d) 2825 active[poly->var] = 1; 2826 2827 rec = isl_poly_as_rec(poly); 2828 for (i = 0; i < rec->n; ++i) 2829 if (poly_set_active(rec->p[i], active, d) < 0) 2830 return isl_stat_error; 2831 2832 return isl_stat_ok; 2833} 2834 2835static isl_stat set_active(__isl_keep isl_qpolynomial *qp, int *active) 2836{ 2837 int i, j; 2838 isl_size d; 2839 isl_space *space; 2840 2841 space = isl_qpolynomial_peek_domain_space(qp); 2842 d = isl_space_dim(space, isl_dim_all); 2843 if (d < 0 || !active) 2844 return isl_stat_error; 2845 2846 for (i = 0; i < d; ++i) 2847 for (j = 0; j < qp->div->n_row; ++j) { 2848 if (isl_int_is_zero(qp->div->row[j][2 + i])) 2849 continue; 2850 active[i] = 1; 2851 break; 2852 } 2853 2854 return poly_set_active(qp->poly, active, d); 2855} 2856 2857#undef TYPE 2858#define TYPE isl_qpolynomial 2859static 2860#include "check_type_range_templ.c" 2861 2862isl_bool isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp, 2863 enum isl_dim_type type, unsigned first, unsigned n) 2864{ 2865 int i; 2866 int *active = NULL; 2867 isl_bool involves = isl_bool_false; 2868 isl_size offset; 2869 isl_size d; 2870 isl_space *space; 2871 2872 if (!qp) 2873 return isl_bool_error; 2874 if (n == 0) 2875 return isl_bool_false; 2876 2877 if (isl_qpolynomial_check_range(qp, type, first, n) < 0) 2878 return isl_bool_error; 2879 isl_assert(qp->dim->ctx, type == isl_dim_param || 2880 type == isl_dim_in, return isl_bool_error); 2881 2882 space = isl_qpolynomial_peek_domain_space(qp); 2883 d = isl_space_dim(space, isl_dim_all); 2884 if (d < 0) 2885 return isl_bool_error; 2886 active = isl_calloc_array(qp->dim->ctx, int, d); 2887 if (set_active(qp, active) < 0) 2888 goto error; 2889 2890 offset = isl_qpolynomial_domain_var_offset(qp, domain_type(type)); 2891 if (offset < 0) 2892 goto error; 2893 first += offset; 2894 for (i = 0; i < n; ++i) 2895 if (active[first + i]) { 2896 involves = isl_bool_true; 2897 break; 2898 } 2899 2900 free(active); 2901 2902 return involves; 2903error: 2904 free(active); 2905 return isl_bool_error; 2906} 2907 2908/* Remove divs that do not appear in the quasi-polynomial, nor in any 2909 * of the divs that do appear in the quasi-polynomial. 2910 */ 2911static __isl_give isl_qpolynomial *remove_redundant_divs( 2912 __isl_take isl_qpolynomial *qp) 2913{ 2914 int i, j; 2915 isl_size div_pos; 2916 int len; 2917 int skip; 2918 int *active = NULL; 2919 int *reordering = NULL; 2920 int redundant = 0; 2921 int n_div; 2922 isl_ctx *ctx; 2923 2924 if (!qp) 2925 return NULL; 2926 if (qp->div->n_row == 0) 2927 return qp; 2928 2929 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 2930 if (div_pos < 0) 2931 return isl_qpolynomial_free(qp); 2932 len = qp->div->n_col - 2; 2933 ctx = isl_qpolynomial_get_ctx(qp); 2934 active = isl_calloc_array(ctx, int, len); 2935 if (!active) 2936 goto error; 2937 2938 if (poly_set_active(qp->poly, active, len) < 0) 2939 goto error; 2940 2941 for (i = qp->div->n_row - 1; i >= 0; --i) { 2942 if (!active[div_pos + i]) { 2943 redundant = 1; 2944 continue; 2945 } 2946 for (j = 0; j < i; ++j) { 2947 if (isl_int_is_zero(qp->div->row[i][2 + div_pos + j])) 2948 continue; 2949 active[div_pos + j] = 1; 2950 break; 2951 } 2952 } 2953 2954 if (!redundant) { 2955 free(active); 2956 return qp; 2957 } 2958 2959 reordering = isl_alloc_array(qp->div->ctx, int, len); 2960 if (!reordering) 2961 goto error; 2962 2963 for (i = 0; i < div_pos; ++i) 2964 reordering[i] = i; 2965 2966 skip = 0; 2967 n_div = qp->div->n_row; 2968 for (i = 0; i < n_div; ++i) { 2969 if (!active[div_pos + i]) { 2970 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1); 2971 qp->div = isl_mat_drop_cols(qp->div, 2972 2 + div_pos + i - skip, 1); 2973 skip++; 2974 } 2975 reordering[div_pos + i] = div_pos + i - skip; 2976 } 2977 2978 qp->poly = reorder(qp->poly, reordering); 2979 2980 if (!qp->poly || !qp->div) 2981 goto error; 2982 2983 free(active); 2984 free(reordering); 2985 2986 return qp; 2987error: 2988 free(active); 2989 free(reordering); 2990 isl_qpolynomial_free(qp); 2991 return NULL; 2992} 2993 2994__isl_give isl_poly *isl_poly_drop(__isl_take isl_poly *poly, 2995 unsigned first, unsigned n) 2996{ 2997 int i; 2998 isl_poly_rec *rec; 2999 3000 if (!poly) 3001 return NULL; 3002 if (n == 0 || poly->var < 0 || poly->var < first) 3003 return poly; 3004 if (poly->var < first + n) { 3005 poly = replace_by_constant_term(poly); 3006 return isl_poly_drop(poly, first, n); 3007 } 3008 poly = isl_poly_cow(poly); 3009 if (!poly) 3010 return NULL; 3011 poly->var -= n; 3012 rec = isl_poly_as_rec(poly); 3013 if (!rec) 3014 goto error; 3015 3016 for (i = 0; i < rec->n; ++i) { 3017 rec->p[i] = isl_poly_drop(rec->p[i], first, n); 3018 if (!rec->p[i]) 3019 goto error; 3020 } 3021 3022 return poly; 3023error: 3024 isl_poly_free(poly); 3025 return NULL; 3026} 3027 3028__isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name( 3029 __isl_take isl_qpolynomial *qp, 3030 enum isl_dim_type type, unsigned pos, const char *s) 3031{ 3032 isl_space *space; 3033 3034 if (!qp) 3035 return NULL; 3036 if (type == isl_dim_out) 3037 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid, 3038 "cannot set name of output/set dimension", 3039 return isl_qpolynomial_free(qp)); 3040 type = domain_type(type); 3041 space = isl_qpolynomial_take_domain_space(qp); 3042 space = isl_space_set_dim_name(space, type, pos, s); 3043 qp = isl_qpolynomial_restore_domain_space(qp, space); 3044 return qp; 3045} 3046 3047__isl_give isl_qpolynomial *isl_qpolynomial_drop_dims( 3048 __isl_take isl_qpolynomial *qp, 3049 enum isl_dim_type type, unsigned first, unsigned n) 3050{ 3051 isl_space *space; 3052 isl_size offset; 3053 3054 if (!qp) 3055 return NULL; 3056 if (type == isl_dim_out) 3057 isl_die(qp->dim->ctx, isl_error_invalid, 3058 "cannot drop output/set dimension", 3059 goto error); 3060 if (isl_qpolynomial_check_range(qp, type, first, n) < 0) 3061 return isl_qpolynomial_free(qp); 3062 type = domain_type(type); 3063 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type)) 3064 return qp; 3065 3066 3067 isl_assert(qp->dim->ctx, type == isl_dim_param || 3068 type == isl_dim_set, goto error); 3069 3070 space = isl_qpolynomial_take_domain_space(qp); 3071 space = isl_space_drop_dims(space, type, first, n); 3072 qp = isl_qpolynomial_restore_domain_space(qp, space); 3073 3074 qp = isl_qpolynomial_cow(qp); 3075 if (!qp) 3076 return NULL; 3077 3078 offset = isl_qpolynomial_domain_var_offset(qp, type); 3079 if (offset < 0) 3080 goto error; 3081 first += offset; 3082 3083 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n); 3084 if (!qp->div) 3085 goto error; 3086 3087 qp->poly = isl_poly_drop(qp->poly, first, n); 3088 if (!qp->poly) 3089 goto error; 3090 3091 return qp; 3092error: 3093 isl_qpolynomial_free(qp); 3094 return NULL; 3095} 3096 3097/* Project the domain of the quasi-polynomial onto its parameter space. 3098 * The quasi-polynomial may not involve any of the domain dimensions. 3099 */ 3100__isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params( 3101 __isl_take isl_qpolynomial *qp) 3102{ 3103 isl_space *space; 3104 isl_size n; 3105 isl_bool involves; 3106 3107 n = isl_qpolynomial_dim(qp, isl_dim_in); 3108 if (n < 0) 3109 return isl_qpolynomial_free(qp); 3110 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n); 3111 if (involves < 0) 3112 return isl_qpolynomial_free(qp); 3113 if (involves) 3114 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid, 3115 "polynomial involves some of the domain dimensions", 3116 return isl_qpolynomial_free(qp)); 3117 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n); 3118 space = isl_qpolynomial_get_domain_space(qp); 3119 space = isl_space_params(space); 3120 qp = isl_qpolynomial_reset_domain_space(qp, space); 3121 return qp; 3122} 3123 3124static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted( 3125 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq) 3126{ 3127 int i, j, k; 3128 isl_int denom; 3129 unsigned total; 3130 unsigned n_div; 3131 isl_poly *poly; 3132 3133 if (!eq) 3134 goto error; 3135 if (eq->n_eq == 0) { 3136 isl_basic_set_free(eq); 3137 return qp; 3138 } 3139 3140 qp = isl_qpolynomial_cow(qp); 3141 if (!qp) 3142 goto error; 3143 qp->div = isl_mat_cow(qp->div); 3144 if (!qp->div) 3145 goto error; 3146 3147 total = isl_basic_set_offset(eq, isl_dim_div); 3148 n_div = eq->n_div; 3149 isl_int_init(denom); 3150 for (i = 0; i < eq->n_eq; ++i) { 3151 j = isl_seq_last_non_zero(eq->eq[i], total + n_div); 3152 if (j < 0 || j == 0 || j >= total) 3153 continue; 3154 3155 for (k = 0; k < qp->div->n_row; ++k) { 3156 if (isl_int_is_zero(qp->div->row[k][1 + j])) 3157 continue; 3158 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total, 3159 &qp->div->row[k][0]); 3160 normalize_div(qp, k); 3161 } 3162 3163 if (isl_int_is_pos(eq->eq[i][j])) 3164 isl_seq_neg(eq->eq[i], eq->eq[i], total); 3165 isl_int_abs(denom, eq->eq[i][j]); 3166 isl_int_set_si(eq->eq[i][j], 0); 3167 3168 poly = isl_poly_from_affine(qp->dim->ctx, 3169 eq->eq[i], denom, total); 3170 qp->poly = isl_poly_subs(qp->poly, j - 1, 1, &poly); 3171 isl_poly_free(poly); 3172 } 3173 isl_int_clear(denom); 3174 3175 if (!qp->poly) 3176 goto error; 3177 3178 isl_basic_set_free(eq); 3179 3180 qp = substitute_non_divs(qp); 3181 qp = sort_divs(qp); 3182 3183 return qp; 3184error: 3185 isl_basic_set_free(eq); 3186 isl_qpolynomial_free(qp); 3187 return NULL; 3188} 3189 3190/* Exploit the equalities in "eq" to simplify the quasi-polynomial. 3191 */ 3192__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities( 3193 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq) 3194{ 3195 if (!qp || !eq) 3196 goto error; 3197 if (qp->div->n_row > 0) 3198 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row); 3199 return isl_qpolynomial_substitute_equalities_lifted(qp, eq); 3200error: 3201 isl_basic_set_free(eq); 3202 isl_qpolynomial_free(qp); 3203 return NULL; 3204} 3205 3206/* Look for equalities among the variables shared by context and qp 3207 * and the integer divisions of qp, if any. 3208 * The equalities are then used to eliminate variables and/or integer 3209 * divisions from qp. 3210 */ 3211__isl_give isl_qpolynomial *isl_qpolynomial_gist( 3212 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context) 3213{ 3214 isl_local_space *ls; 3215 isl_basic_set *aff; 3216 3217 ls = isl_qpolynomial_get_domain_local_space(qp); 3218 context = isl_local_space_lift_set(ls, context); 3219 3220 aff = isl_set_affine_hull(context); 3221 return isl_qpolynomial_substitute_equalities_lifted(qp, aff); 3222} 3223 3224__isl_give isl_qpolynomial *isl_qpolynomial_gist_params( 3225 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context) 3226{ 3227 isl_space *space = isl_qpolynomial_get_domain_space(qp); 3228 isl_set *dom_context = isl_set_universe(space); 3229 dom_context = isl_set_intersect_params(dom_context, context); 3230 return isl_qpolynomial_gist(qp, dom_context); 3231} 3232 3233/* Return a zero isl_qpolynomial in the given space. 3234 * 3235 * This is a helper function for isl_pw_*_as_* that ensures a uniform 3236 * interface over all piecewise types. 3237 */ 3238static __isl_give isl_qpolynomial *isl_qpolynomial_zero_in_space( 3239 __isl_take isl_space *space) 3240{ 3241 return isl_qpolynomial_zero_on_domain(isl_space_domain(space)); 3242} 3243 3244#define isl_qpolynomial_involves_nan isl_qpolynomial_is_nan 3245 3246#undef PW 3247#define PW isl_pw_qpolynomial 3248#undef BASE 3249#define BASE qpolynomial 3250#undef EL_IS_ZERO 3251#define EL_IS_ZERO is_zero 3252#undef ZERO 3253#define ZERO zero 3254#undef IS_ZERO 3255#define IS_ZERO is_zero 3256#undef FIELD 3257#define FIELD qp 3258#undef DEFAULT_IS_ZERO 3259#define DEFAULT_IS_ZERO 1 3260 3261#include <isl_pw_templ.c> 3262#include <isl_pw_un_op_templ.c> 3263#include <isl_pw_add_disjoint_templ.c> 3264#include <isl_pw_domain_reverse_templ.c> 3265#include <isl_pw_eval.c> 3266#include <isl_pw_fix_templ.c> 3267#include <isl_pw_from_range_templ.c> 3268#include <isl_pw_insert_dims_templ.c> 3269#include <isl_pw_lift_templ.c> 3270#include <isl_pw_morph_templ.c> 3271#include <isl_pw_move_dims_templ.c> 3272#include <isl_pw_neg_templ.c> 3273#include <isl_pw_opt_templ.c> 3274#include <isl_pw_split_dims_templ.c> 3275#include <isl_pw_sub_templ.c> 3276 3277#undef BASE 3278#define BASE pw_qpolynomial 3279 3280#include <isl_union_single.c> 3281#include <isl_union_domain_reverse_templ.c> 3282#include <isl_union_eval.c> 3283#include <isl_union_neg.c> 3284#include <isl_union_sub_templ.c> 3285 3286int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp) 3287{ 3288 if (!pwqp) 3289 return -1; 3290 3291 if (pwqp->n != -1) 3292 return 0; 3293 3294 if (!isl_set_plain_is_universe(pwqp->p[0].set)) 3295 return 0; 3296 3297 return isl_qpolynomial_is_one(pwqp->p[0].qp); 3298} 3299 3300__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add( 3301 __isl_take isl_pw_qpolynomial *pwqp1, 3302 __isl_take isl_pw_qpolynomial *pwqp2) 3303{ 3304 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2); 3305} 3306 3307__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul( 3308 __isl_take isl_pw_qpolynomial *pwqp1, 3309 __isl_take isl_pw_qpolynomial *pwqp2) 3310{ 3311 int i, j, n; 3312 struct isl_pw_qpolynomial *res; 3313 3314 if (!pwqp1 || !pwqp2) 3315 goto error; 3316 3317 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim), 3318 goto error); 3319 3320 if (isl_pw_qpolynomial_is_zero(pwqp1)) { 3321 isl_pw_qpolynomial_free(pwqp2); 3322 return pwqp1; 3323 } 3324 3325 if (isl_pw_qpolynomial_is_zero(pwqp2)) { 3326 isl_pw_qpolynomial_free(pwqp1); 3327 return pwqp2; 3328 } 3329 3330 if (isl_pw_qpolynomial_is_one(pwqp1)) { 3331 isl_pw_qpolynomial_free(pwqp1); 3332 return pwqp2; 3333 } 3334 3335 if (isl_pw_qpolynomial_is_one(pwqp2)) { 3336 isl_pw_qpolynomial_free(pwqp2); 3337 return pwqp1; 3338 } 3339 3340 n = pwqp1->n * pwqp2->n; 3341 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n); 3342 3343 for (i = 0; i < pwqp1->n; ++i) { 3344 for (j = 0; j < pwqp2->n; ++j) { 3345 struct isl_set *common; 3346 struct isl_qpolynomial *prod; 3347 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set), 3348 isl_set_copy(pwqp2->p[j].set)); 3349 if (isl_set_plain_is_empty(common)) { 3350 isl_set_free(common); 3351 continue; 3352 } 3353 3354 prod = isl_qpolynomial_mul( 3355 isl_qpolynomial_copy(pwqp1->p[i].qp), 3356 isl_qpolynomial_copy(pwqp2->p[j].qp)); 3357 3358 res = isl_pw_qpolynomial_add_piece(res, common, prod); 3359 } 3360 } 3361 3362 isl_pw_qpolynomial_free(pwqp1); 3363 isl_pw_qpolynomial_free(pwqp2); 3364 3365 return res; 3366error: 3367 isl_pw_qpolynomial_free(pwqp1); 3368 isl_pw_qpolynomial_free(pwqp2); 3369 return NULL; 3370} 3371 3372__isl_give isl_val *isl_poly_eval(__isl_take isl_poly *poly, 3373 __isl_take isl_vec *vec) 3374{ 3375 int i; 3376 isl_bool is_cst; 3377 isl_poly_rec *rec; 3378 isl_val *res; 3379 isl_val *base; 3380 3381 is_cst = isl_poly_is_cst(poly); 3382 if (is_cst < 0) 3383 goto error; 3384 if (is_cst) { 3385 isl_vec_free(vec); 3386 res = isl_poly_get_constant_val(poly); 3387 isl_poly_free(poly); 3388 return res; 3389 } 3390 3391 rec = isl_poly_as_rec(poly); 3392 if (!rec || !vec) 3393 goto error; 3394 3395 isl_assert(poly->ctx, rec->n >= 1, goto error); 3396 3397 base = isl_val_rat_from_isl_int(poly->ctx, 3398 vec->el[1 + poly->var], vec->el[0]); 3399 3400 res = isl_poly_eval(isl_poly_copy(rec->p[rec->n - 1]), 3401 isl_vec_copy(vec)); 3402 3403 for (i = rec->n - 2; i >= 0; --i) { 3404 res = isl_val_mul(res, isl_val_copy(base)); 3405 res = isl_val_add(res, isl_poly_eval(isl_poly_copy(rec->p[i]), 3406 isl_vec_copy(vec))); 3407 } 3408 3409 isl_val_free(base); 3410 isl_poly_free(poly); 3411 isl_vec_free(vec); 3412 return res; 3413error: 3414 isl_poly_free(poly); 3415 isl_vec_free(vec); 3416 return NULL; 3417} 3418 3419/* Evaluate "qp" in the void point "pnt". 3420 * In particular, return the value NaN. 3421 */ 3422static __isl_give isl_val *eval_void(__isl_take isl_qpolynomial *qp, 3423 __isl_take isl_point *pnt) 3424{ 3425 isl_ctx *ctx; 3426 3427 ctx = isl_point_get_ctx(pnt); 3428 isl_qpolynomial_free(qp); 3429 isl_point_free(pnt); 3430 return isl_val_nan(ctx); 3431} 3432 3433__isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp, 3434 __isl_take isl_point *pnt) 3435{ 3436 isl_bool is_void; 3437 isl_vec *ext; 3438 isl_val *v; 3439 3440 if (!qp || !pnt) 3441 goto error; 3442 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error); 3443 is_void = isl_point_is_void(pnt); 3444 if (is_void < 0) 3445 goto error; 3446 if (is_void) 3447 return eval_void(qp, pnt); 3448 3449 ext = isl_local_extend_point_vec(qp->div, isl_vec_copy(pnt->vec)); 3450 3451 v = isl_poly_eval(isl_qpolynomial_get_poly(qp), ext); 3452 3453 isl_qpolynomial_free(qp); 3454 isl_point_free(pnt); 3455 3456 return v; 3457error: 3458 isl_qpolynomial_free(qp); 3459 isl_point_free(pnt); 3460 return NULL; 3461} 3462 3463int isl_poly_cmp(__isl_keep isl_poly_cst *cst1, __isl_keep isl_poly_cst *cst2) 3464{ 3465 int cmp; 3466 isl_int t; 3467 isl_int_init(t); 3468 isl_int_mul(t, cst1->n, cst2->d); 3469 isl_int_submul(t, cst2->n, cst1->d); 3470 cmp = isl_int_sgn(t); 3471 isl_int_clear(t); 3472 return cmp; 3473} 3474 3475__isl_give isl_qpolynomial *isl_qpolynomial_insert_dims( 3476 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, 3477 unsigned first, unsigned n) 3478{ 3479 unsigned total; 3480 unsigned g_pos; 3481 int *exp; 3482 isl_space *space; 3483 3484 if (!qp) 3485 return NULL; 3486 if (type == isl_dim_out) 3487 isl_die(qp->div->ctx, isl_error_invalid, 3488 "cannot insert output/set dimensions", 3489 goto error); 3490 if (isl_qpolynomial_check_range(qp, type, first, 0) < 0) 3491 return isl_qpolynomial_free(qp); 3492 type = domain_type(type); 3493 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type)) 3494 return qp; 3495 3496 qp = isl_qpolynomial_cow(qp); 3497 if (!qp) 3498 return NULL; 3499 3500 g_pos = pos(qp->dim, type) + first; 3501 3502 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n); 3503 if (!qp->div) 3504 goto error; 3505 3506 total = qp->div->n_col - 2; 3507 if (total > g_pos) { 3508 int i; 3509 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos); 3510 if (!exp) 3511 goto error; 3512 for (i = 0; i < total - g_pos; ++i) 3513 exp[i] = i + n; 3514 qp->poly = expand(qp->poly, exp, g_pos); 3515 free(exp); 3516 if (!qp->poly) 3517 goto error; 3518 } 3519 3520 space = isl_qpolynomial_take_domain_space(qp); 3521 space = isl_space_insert_dims(space, type, first, n); 3522 qp = isl_qpolynomial_restore_domain_space(qp, space); 3523 3524 return qp; 3525error: 3526 isl_qpolynomial_free(qp); 3527 return NULL; 3528} 3529 3530__isl_give isl_qpolynomial *isl_qpolynomial_add_dims( 3531 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n) 3532{ 3533 isl_size pos; 3534 3535 pos = isl_qpolynomial_dim(qp, type); 3536 if (pos < 0) 3537 return isl_qpolynomial_free(qp); 3538 3539 return isl_qpolynomial_insert_dims(qp, type, pos, n); 3540} 3541 3542static int *reordering_move(isl_ctx *ctx, 3543 unsigned len, unsigned dst, unsigned src, unsigned n) 3544{ 3545 int i; 3546 int *reordering; 3547 3548 reordering = isl_alloc_array(ctx, int, len); 3549 if (!reordering) 3550 return NULL; 3551 3552 if (dst <= src) { 3553 for (i = 0; i < dst; ++i) 3554 reordering[i] = i; 3555 for (i = 0; i < n; ++i) 3556 reordering[src + i] = dst + i; 3557 for (i = 0; i < src - dst; ++i) 3558 reordering[dst + i] = dst + n + i; 3559 for (i = 0; i < len - src - n; ++i) 3560 reordering[src + n + i] = src + n + i; 3561 } else { 3562 for (i = 0; i < src; ++i) 3563 reordering[i] = i; 3564 for (i = 0; i < n; ++i) 3565 reordering[src + i] = dst + i; 3566 for (i = 0; i < dst - src; ++i) 3567 reordering[src + n + i] = src + i; 3568 for (i = 0; i < len - dst - n; ++i) 3569 reordering[dst + n + i] = dst + n + i; 3570 } 3571 3572 return reordering; 3573} 3574 3575/* Move the "n" variables starting at "src_pos" of "qp" to "dst_pos". 3576 * Only modify the polynomial expression and the local variables of "qp". 3577 * The caller is responsible for modifying the space accordingly. 3578 */ 3579static __isl_give isl_qpolynomial *local_poly_move_dims( 3580 __isl_take isl_qpolynomial *qp, 3581 unsigned dst_pos, unsigned src_pos, unsigned n) 3582{ 3583 isl_ctx *ctx; 3584 isl_size total; 3585 int *reordering; 3586 isl_local *local; 3587 isl_poly *poly; 3588 3589 local = isl_qpolynomial_take_local(qp); 3590 local = isl_local_move_vars(local, dst_pos, src_pos, n); 3591 qp = isl_qpolynomial_restore_local(qp, local); 3592 qp = sort_divs(qp); 3593 3594 total = isl_qpolynomial_domain_dim(qp, isl_dim_all); 3595 if (total < 0) 3596 return isl_qpolynomial_free(qp); 3597 ctx = isl_qpolynomial_get_ctx(qp); 3598 reordering = reordering_move(ctx, total, dst_pos, src_pos, n); 3599 if (!reordering) 3600 return isl_qpolynomial_free(qp); 3601 3602 poly = isl_qpolynomial_take_poly(qp); 3603 poly = reorder(poly, reordering); 3604 qp = isl_qpolynomial_restore_poly(qp, poly); 3605 free(reordering); 3606 3607 return qp; 3608} 3609 3610__isl_give isl_qpolynomial *isl_qpolynomial_move_dims( 3611 __isl_take isl_qpolynomial *qp, 3612 enum isl_dim_type dst_type, unsigned dst_pos, 3613 enum isl_dim_type src_type, unsigned src_pos, unsigned n) 3614{ 3615 isl_ctx *ctx; 3616 unsigned g_dst_pos; 3617 unsigned g_src_pos; 3618 isl_size src_off, dst_off; 3619 isl_space *space; 3620 3621 if (!qp) 3622 return NULL; 3623 3624 ctx = isl_qpolynomial_get_ctx(qp); 3625 if (dst_type == isl_dim_out || src_type == isl_dim_out) 3626 isl_die(ctx, isl_error_invalid, 3627 "cannot move output/set dimension", 3628 return isl_qpolynomial_free(qp)); 3629 if (src_type == isl_dim_div || dst_type == isl_dim_div) 3630 isl_die(ctx, isl_error_invalid, "cannot move local variables", 3631 return isl_qpolynomial_free(qp)); 3632 if (isl_qpolynomial_check_range(qp, src_type, src_pos, n) < 0) 3633 return isl_qpolynomial_free(qp); 3634 if (dst_type == isl_dim_in) 3635 dst_type = isl_dim_set; 3636 if (src_type == isl_dim_in) 3637 src_type = isl_dim_set; 3638 3639 if (n == 0 && 3640 !isl_space_is_named_or_nested(qp->dim, src_type) && 3641 !isl_space_is_named_or_nested(qp->dim, dst_type)) 3642 return qp; 3643 3644 src_off = isl_qpolynomial_domain_var_offset(qp, src_type); 3645 dst_off = isl_qpolynomial_domain_var_offset(qp, dst_type); 3646 if (src_off < 0 || dst_off < 0) 3647 return isl_qpolynomial_free(qp); 3648 3649 g_dst_pos = dst_off + dst_pos; 3650 g_src_pos = src_off + src_pos; 3651 if (dst_type > src_type) 3652 g_dst_pos -= n; 3653 3654 qp = local_poly_move_dims(qp, g_dst_pos, g_src_pos, n); 3655 3656 space = isl_qpolynomial_take_domain_space(qp); 3657 space = isl_space_move_dims(space, dst_type, dst_pos, 3658 src_type, src_pos, n); 3659 qp = isl_qpolynomial_restore_domain_space(qp, space); 3660 3661 return qp; 3662} 3663 3664/* Given a quasi-polynomial on a domain (A -> B), 3665 * interchange A and B in the wrapped domain 3666 * to obtain a quasi-polynomial on the domain (B -> A). 3667 */ 3668__isl_give isl_qpolynomial *isl_qpolynomial_domain_reverse( 3669 __isl_take isl_qpolynomial *qp) 3670{ 3671 isl_space *space; 3672 isl_size n_in, n_out, offset; 3673 3674 space = isl_qpolynomial_peek_domain_space(qp); 3675 offset = isl_space_offset(space, isl_dim_set); 3676 n_in = isl_space_wrapped_dim(space, isl_dim_set, isl_dim_in); 3677 n_out = isl_space_wrapped_dim(space, isl_dim_set, isl_dim_out); 3678 if (offset < 0 || n_in < 0 || n_out < 0) 3679 return isl_qpolynomial_free(qp); 3680 3681 qp = local_poly_move_dims(qp, offset, offset + n_in, n_out); 3682 3683 space = isl_qpolynomial_take_domain_space(qp); 3684 space = isl_space_wrapped_reverse(space); 3685 qp = isl_qpolynomial_restore_domain_space(qp, space); 3686 3687 return qp; 3688} 3689 3690__isl_give isl_qpolynomial *isl_qpolynomial_from_affine( 3691 __isl_take isl_space *space, isl_int *f, isl_int denom) 3692{ 3693 isl_size d; 3694 isl_poly *poly; 3695 3696 space = isl_space_domain(space); 3697 if (!space) 3698 return NULL; 3699 3700 d = isl_space_dim(space, isl_dim_all); 3701 poly = d < 0 ? NULL : isl_poly_from_affine(space->ctx, f, denom, 1 + d); 3702 3703 return isl_qpolynomial_alloc(space, 0, poly); 3704} 3705 3706__isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff) 3707{ 3708 isl_ctx *ctx; 3709 isl_poly *poly; 3710 isl_qpolynomial *qp; 3711 3712 if (!aff) 3713 return NULL; 3714 3715 ctx = isl_aff_get_ctx(aff); 3716 poly = isl_poly_from_affine(ctx, aff->v->el + 1, aff->v->el[0], 3717 aff->v->size - 1); 3718 3719 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff), 3720 aff->ls->div->n_row, poly); 3721 if (!qp) 3722 goto error; 3723 3724 isl_mat_free(qp->div); 3725 qp->div = isl_mat_copy(aff->ls->div); 3726 qp->div = isl_mat_cow(qp->div); 3727 if (!qp->div) 3728 goto error; 3729 3730 isl_aff_free(aff); 3731 qp = reduce_divs(qp); 3732 qp = remove_redundant_divs(qp); 3733 return qp; 3734error: 3735 isl_aff_free(aff); 3736 return isl_qpolynomial_free(qp); 3737} 3738 3739__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff( 3740 __isl_take isl_pw_aff *pwaff) 3741{ 3742 int i; 3743 isl_pw_qpolynomial *pwqp; 3744 3745 if (!pwaff) 3746 return NULL; 3747 3748 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff), 3749 pwaff->n); 3750 3751 for (i = 0; i < pwaff->n; ++i) { 3752 isl_set *dom; 3753 isl_qpolynomial *qp; 3754 3755 dom = isl_set_copy(pwaff->p[i].set); 3756 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff)); 3757 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp); 3758 } 3759 3760 isl_pw_aff_free(pwaff); 3761 return pwqp; 3762} 3763 3764__isl_give isl_qpolynomial *isl_qpolynomial_from_constraint( 3765 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos) 3766{ 3767 isl_aff *aff; 3768 3769 aff = isl_constraint_get_bound(c, type, pos); 3770 isl_constraint_free(c); 3771 return isl_qpolynomial_from_aff(aff); 3772} 3773 3774/* For each 0 <= i < "n", replace variable "first" + i of type "type" 3775 * in "qp" by subs[i]. 3776 */ 3777__isl_give isl_qpolynomial *isl_qpolynomial_substitute( 3778 __isl_take isl_qpolynomial *qp, 3779 enum isl_dim_type type, unsigned first, unsigned n, 3780 __isl_keep isl_qpolynomial **subs) 3781{ 3782 int i; 3783 isl_poly *poly; 3784 isl_poly **polys; 3785 3786 if (n == 0) 3787 return qp; 3788 3789 if (!qp) 3790 return NULL; 3791 3792 if (type == isl_dim_out) 3793 isl_die(qp->dim->ctx, isl_error_invalid, 3794 "cannot substitute output/set dimension", 3795 goto error); 3796 if (isl_qpolynomial_check_range(qp, type, first, n) < 0) 3797 return isl_qpolynomial_free(qp); 3798 type = domain_type(type); 3799 3800 for (i = 0; i < n; ++i) 3801 if (!subs[i]) 3802 goto error; 3803 3804 for (i = 0; i < n; ++i) 3805 if (isl_qpolynomial_check_equal_space(qp, subs[i]) < 0) 3806 goto error; 3807 3808 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error); 3809 for (i = 0; i < n; ++i) 3810 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error); 3811 3812 first += pos(qp->dim, type); 3813 3814 polys = isl_alloc_array(qp->dim->ctx, struct isl_poly *, n); 3815 if (!polys) 3816 goto error; 3817 for (i = 0; i < n; ++i) 3818 polys[i] = subs[i]->poly; 3819 3820 poly = isl_qpolynomial_take_poly(qp); 3821 poly = isl_poly_subs(poly, first, n, polys); 3822 qp = isl_qpolynomial_restore_poly(qp, poly); 3823 3824 free(polys); 3825 3826 return qp; 3827error: 3828 isl_qpolynomial_free(qp); 3829 return NULL; 3830} 3831 3832/* Extend "bset" with extra set dimensions for each integer division 3833 * in "qp" and then call "fn" with the extended bset and the polynomial 3834 * that results from replacing each of the integer divisions by the 3835 * corresponding extra set dimension. 3836 */ 3837isl_stat isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp, 3838 __isl_keep isl_basic_set *bset, 3839 isl_stat (*fn)(__isl_take isl_basic_set *bset, 3840 __isl_take isl_qpolynomial *poly, void *user), void *user) 3841{ 3842 isl_space *space; 3843 isl_local_space *ls; 3844 isl_poly *poly; 3845 isl_qpolynomial *polynomial; 3846 3847 if (!qp || !bset) 3848 return isl_stat_error; 3849 if (qp->div->n_row == 0) 3850 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp), 3851 user); 3852 3853 space = isl_space_copy(qp->dim); 3854 space = isl_space_add_dims(space, isl_dim_set, qp->div->n_row); 3855 poly = isl_qpolynomial_get_poly(qp); 3856 polynomial = isl_qpolynomial_alloc(space, 0, poly); 3857 bset = isl_basic_set_copy(bset); 3858 ls = isl_qpolynomial_get_domain_local_space(qp); 3859 bset = isl_local_space_lift_basic_set(ls, bset); 3860 3861 return fn(bset, polynomial, user); 3862} 3863 3864/* Return total degree in variables first (inclusive) up to last (exclusive). 3865 */ 3866int isl_poly_degree(__isl_keep isl_poly *poly, int first, int last) 3867{ 3868 int deg = -1; 3869 int i; 3870 isl_bool is_zero, is_cst; 3871 isl_poly_rec *rec; 3872 3873 is_zero = isl_poly_is_zero(poly); 3874 if (is_zero < 0) 3875 return -2; 3876 if (is_zero) 3877 return -1; 3878 is_cst = isl_poly_is_cst(poly); 3879 if (is_cst < 0) 3880 return -2; 3881 if (is_cst || poly->var < first) 3882 return 0; 3883 3884 rec = isl_poly_as_rec(poly); 3885 if (!rec) 3886 return -2; 3887 3888 for (i = 0; i < rec->n; ++i) { 3889 int d; 3890 3891 is_zero = isl_poly_is_zero(rec->p[i]); 3892 if (is_zero < 0) 3893 return -2; 3894 if (is_zero) 3895 continue; 3896 d = isl_poly_degree(rec->p[i], first, last); 3897 if (poly->var < last) 3898 d += i; 3899 if (d > deg) 3900 deg = d; 3901 } 3902 3903 return deg; 3904} 3905 3906/* Return total degree in set variables. 3907 */ 3908int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly) 3909{ 3910 isl_size ovar; 3911 isl_size nvar; 3912 3913 if (!poly) 3914 return -2; 3915 3916 ovar = isl_space_offset(poly->dim, isl_dim_set); 3917 nvar = isl_space_dim(poly->dim, isl_dim_set); 3918 if (ovar < 0 || nvar < 0) 3919 return -2; 3920 return isl_poly_degree(poly->poly, ovar, ovar + nvar); 3921} 3922 3923__isl_give isl_poly *isl_poly_coeff(__isl_keep isl_poly *poly, 3924 unsigned pos, int deg) 3925{ 3926 int i; 3927 isl_bool is_cst; 3928 isl_poly_rec *rec; 3929 3930 is_cst = isl_poly_is_cst(poly); 3931 if (is_cst < 0) 3932 return NULL; 3933 if (is_cst || poly->var < pos) { 3934 if (deg == 0) 3935 return isl_poly_copy(poly); 3936 else 3937 return isl_poly_zero(poly->ctx); 3938 } 3939 3940 rec = isl_poly_as_rec(poly); 3941 if (!rec) 3942 return NULL; 3943 3944 if (poly->var == pos) { 3945 if (deg < rec->n) 3946 return isl_poly_copy(rec->p[deg]); 3947 else 3948 return isl_poly_zero(poly->ctx); 3949 } 3950 3951 poly = isl_poly_copy(poly); 3952 poly = isl_poly_cow(poly); 3953 rec = isl_poly_as_rec(poly); 3954 if (!rec) 3955 goto error; 3956 3957 for (i = 0; i < rec->n; ++i) { 3958 isl_poly *t; 3959 t = isl_poly_coeff(rec->p[i], pos, deg); 3960 if (!t) 3961 goto error; 3962 isl_poly_free(rec->p[i]); 3963 rec->p[i] = t; 3964 } 3965 3966 return poly; 3967error: 3968 isl_poly_free(poly); 3969 return NULL; 3970} 3971 3972/* Return coefficient of power "deg" of variable "t_pos" of type "type". 3973 */ 3974__isl_give isl_qpolynomial *isl_qpolynomial_coeff( 3975 __isl_keep isl_qpolynomial *qp, 3976 enum isl_dim_type type, unsigned t_pos, int deg) 3977{ 3978 unsigned g_pos; 3979 isl_poly *poly; 3980 isl_qpolynomial *c; 3981 3982 if (!qp) 3983 return NULL; 3984 3985 if (type == isl_dim_out) 3986 isl_die(qp->div->ctx, isl_error_invalid, 3987 "output/set dimension does not have a coefficient", 3988 return NULL); 3989 if (isl_qpolynomial_check_range(qp, type, t_pos, 1) < 0) 3990 return NULL; 3991 type = domain_type(type); 3992 3993 g_pos = pos(qp->dim, type) + t_pos; 3994 poly = isl_poly_coeff(qp->poly, g_pos, deg); 3995 3996 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), 3997 qp->div->n_row, poly); 3998 if (!c) 3999 return NULL; 4000 isl_mat_free(c->div); 4001 c->div = isl_qpolynomial_get_local(qp); 4002 if (!c->div) 4003 goto error; 4004 return c; 4005error: 4006 isl_qpolynomial_free(c); 4007 return NULL; 4008} 4009 4010/* Homogenize the polynomial in the variables first (inclusive) up to 4011 * last (exclusive) by inserting powers of variable first. 4012 * Variable first is assumed not to appear in the input. 4013 */ 4014__isl_give isl_poly *isl_poly_homogenize(__isl_take isl_poly *poly, int deg, 4015 int target, int first, int last) 4016{ 4017 int i; 4018 isl_bool is_zero, is_cst; 4019 isl_poly_rec *rec; 4020 4021 is_zero = isl_poly_is_zero(poly); 4022 if (is_zero < 0) 4023 return isl_poly_free(poly); 4024 if (is_zero) 4025 return poly; 4026 if (deg == target) 4027 return poly; 4028 is_cst = isl_poly_is_cst(poly); 4029 if (is_cst < 0) 4030 return isl_poly_free(poly); 4031 if (is_cst || poly->var < first) { 4032 isl_poly *hom; 4033 4034 hom = isl_poly_var_pow(poly->ctx, first, target - deg); 4035 if (!hom) 4036 goto error; 4037 rec = isl_poly_as_rec(hom); 4038 rec->p[target - deg] = isl_poly_mul(rec->p[target - deg], poly); 4039 4040 return hom; 4041 } 4042 4043 poly = isl_poly_cow(poly); 4044 rec = isl_poly_as_rec(poly); 4045 if (!rec) 4046 goto error; 4047 4048 for (i = 0; i < rec->n; ++i) { 4049 is_zero = isl_poly_is_zero(rec->p[i]); 4050 if (is_zero < 0) 4051 return isl_poly_free(poly); 4052 if (is_zero) 4053 continue; 4054 rec->p[i] = isl_poly_homogenize(rec->p[i], 4055 poly->var < last ? deg + i : i, target, 4056 first, last); 4057 if (!rec->p[i]) 4058 goto error; 4059 } 4060 4061 return poly; 4062error: 4063 isl_poly_free(poly); 4064 return NULL; 4065} 4066 4067/* Homogenize the polynomial in the set variables by introducing 4068 * powers of an extra set variable at position 0. 4069 */ 4070__isl_give isl_qpolynomial *isl_qpolynomial_homogenize( 4071 __isl_take isl_qpolynomial *poly) 4072{ 4073 isl_size ovar; 4074 isl_size nvar; 4075 int deg = isl_qpolynomial_degree(poly); 4076 4077 if (deg < -1) 4078 goto error; 4079 4080 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1); 4081 poly = isl_qpolynomial_cow(poly); 4082 if (!poly) 4083 goto error; 4084 4085 ovar = isl_space_offset(poly->dim, isl_dim_set); 4086 nvar = isl_space_dim(poly->dim, isl_dim_set); 4087 if (ovar < 0 || nvar < 0) 4088 return isl_qpolynomial_free(poly); 4089 poly->poly = isl_poly_homogenize(poly->poly, 0, deg, ovar, ovar + nvar); 4090 if (!poly->poly) 4091 goto error; 4092 4093 return poly; 4094error: 4095 isl_qpolynomial_free(poly); 4096 return NULL; 4097} 4098 4099__isl_give isl_term *isl_term_alloc(__isl_take isl_space *space, 4100 __isl_take isl_mat *div) 4101{ 4102 isl_term *term; 4103 isl_size d; 4104 int n; 4105 4106 d = isl_space_dim(space, isl_dim_all); 4107 if (d < 0 || !div) 4108 goto error; 4109 4110 n = d + div->n_row; 4111 4112 term = isl_calloc(space->ctx, struct isl_term, 4113 sizeof(struct isl_term) + (n - 1) * sizeof(int)); 4114 if (!term) 4115 goto error; 4116 4117 term->ref = 1; 4118 term->dim = space; 4119 term->div = div; 4120 isl_int_init(term->n); 4121 isl_int_init(term->d); 4122 4123 return term; 4124error: 4125 isl_space_free(space); 4126 isl_mat_free(div); 4127 return NULL; 4128} 4129 4130__isl_give isl_term *isl_term_copy(__isl_keep isl_term *term) 4131{ 4132 if (!term) 4133 return NULL; 4134 4135 term->ref++; 4136 return term; 4137} 4138 4139__isl_give isl_term *isl_term_dup(__isl_keep isl_term *term) 4140{ 4141 int i; 4142 isl_term *dup; 4143 isl_size total; 4144 4145 total = isl_term_dim(term, isl_dim_all); 4146 if (total < 0) 4147 return NULL; 4148 4149 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div)); 4150 if (!dup) 4151 return NULL; 4152 4153 isl_int_set(dup->n, term->n); 4154 isl_int_set(dup->d, term->d); 4155 4156 for (i = 0; i < total; ++i) 4157 dup->pow[i] = term->pow[i]; 4158 4159 return dup; 4160} 4161 4162__isl_give isl_term *isl_term_cow(__isl_take isl_term *term) 4163{ 4164 if (!term) 4165 return NULL; 4166 4167 if (term->ref == 1) 4168 return term; 4169 term->ref--; 4170 return isl_term_dup(term); 4171} 4172 4173__isl_null isl_term *isl_term_free(__isl_take isl_term *term) 4174{ 4175 if (!term) 4176 return NULL; 4177 4178 if (--term->ref > 0) 4179 return NULL; 4180 4181 isl_space_free(term->dim); 4182 isl_mat_free(term->div); 4183 isl_int_clear(term->n); 4184 isl_int_clear(term->d); 4185 free(term); 4186 4187 return NULL; 4188} 4189 4190isl_size isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type) 4191{ 4192 isl_size dim; 4193 4194 if (!term) 4195 return isl_size_error; 4196 4197 switch (type) { 4198 case isl_dim_param: 4199 case isl_dim_in: 4200 case isl_dim_out: return isl_space_dim(term->dim, type); 4201 case isl_dim_div: return term->div->n_row; 4202 case isl_dim_all: dim = isl_space_dim(term->dim, isl_dim_all); 4203 if (dim < 0) 4204 return isl_size_error; 4205 return dim + term->div->n_row; 4206 default: return isl_size_error; 4207 } 4208} 4209 4210/* Return the space of "term". 4211 */ 4212static __isl_keep isl_space *isl_term_peek_space(__isl_keep isl_term *term) 4213{ 4214 return term ? term->dim : NULL; 4215} 4216 4217/* Return the offset of the first variable of type "type" within 4218 * the variables of "term". 4219 */ 4220static isl_size isl_term_offset(__isl_keep isl_term *term, 4221 enum isl_dim_type type) 4222{ 4223 isl_space *space; 4224 4225 space = isl_term_peek_space(term); 4226 if (!space) 4227 return isl_size_error; 4228 4229 switch (type) { 4230 case isl_dim_param: 4231 case isl_dim_set: return isl_space_offset(space, type); 4232 case isl_dim_div: return isl_space_dim(space, isl_dim_all); 4233 default: 4234 isl_die(isl_term_get_ctx(term), isl_error_invalid, 4235 "invalid dimension type", return isl_size_error); 4236 } 4237} 4238 4239isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term) 4240{ 4241 return term ? term->dim->ctx : NULL; 4242} 4243 4244void isl_term_get_num(__isl_keep isl_term *term, isl_int *n) 4245{ 4246 if (!term) 4247 return; 4248 isl_int_set(*n, term->n); 4249} 4250 4251/* Return the coefficient of the term "term". 4252 */ 4253__isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term) 4254{ 4255 if (!term) 4256 return NULL; 4257 4258 return isl_val_rat_from_isl_int(isl_term_get_ctx(term), 4259 term->n, term->d); 4260} 4261 4262#undef TYPE 4263#define TYPE isl_term 4264static 4265#include "check_type_range_templ.c" 4266 4267isl_size isl_term_get_exp(__isl_keep isl_term *term, 4268 enum isl_dim_type type, unsigned pos) 4269{ 4270 isl_size offset; 4271 4272 if (isl_term_check_range(term, type, pos, 1) < 0) 4273 return isl_size_error; 4274 offset = isl_term_offset(term, type); 4275 if (offset < 0) 4276 return isl_size_error; 4277 4278 return term->pow[offset + pos]; 4279} 4280 4281__isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos) 4282{ 4283 isl_local_space *ls; 4284 isl_aff *aff; 4285 4286 if (isl_term_check_range(term, isl_dim_div, pos, 1) < 0) 4287 return NULL; 4288 4289 ls = isl_local_space_alloc_div(isl_space_copy(term->dim), 4290 isl_mat_copy(term->div)); 4291 aff = isl_aff_alloc(ls); 4292 if (!aff) 4293 return NULL; 4294 4295 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size); 4296 4297 aff = isl_aff_normalize(aff); 4298 4299 return aff; 4300} 4301 4302__isl_give isl_term *isl_poly_foreach_term(__isl_keep isl_poly *poly, 4303 isl_stat (*fn)(__isl_take isl_term *term, void *user), 4304 __isl_take isl_term *term, void *user) 4305{ 4306 int i; 4307 isl_bool is_zero, is_bad, is_cst; 4308 isl_poly_rec *rec; 4309 4310 is_zero = isl_poly_is_zero(poly); 4311 if (is_zero < 0 || !term) 4312 goto error; 4313 4314 if (is_zero) 4315 return term; 4316 4317 is_cst = isl_poly_is_cst(poly); 4318 is_bad = isl_poly_is_nan(poly); 4319 if (is_bad >= 0 && !is_bad) 4320 is_bad = isl_poly_is_infty(poly); 4321 if (is_bad >= 0 && !is_bad) 4322 is_bad = isl_poly_is_neginfty(poly); 4323 if (is_cst < 0 || is_bad < 0) 4324 return isl_term_free(term); 4325 if (is_bad) 4326 isl_die(isl_term_get_ctx(term), isl_error_invalid, 4327 "cannot handle NaN/infty polynomial", 4328 return isl_term_free(term)); 4329 4330 if (is_cst) { 4331 isl_poly_cst *cst; 4332 cst = isl_poly_as_cst(poly); 4333 if (!cst) 4334 goto error; 4335 term = isl_term_cow(term); 4336 if (!term) 4337 goto error; 4338 isl_int_set(term->n, cst->n); 4339 isl_int_set(term->d, cst->d); 4340 if (fn(isl_term_copy(term), user) < 0) 4341 goto error; 4342 return term; 4343 } 4344 4345 rec = isl_poly_as_rec(poly); 4346 if (!rec) 4347 goto error; 4348 4349 for (i = 0; i < rec->n; ++i) { 4350 term = isl_term_cow(term); 4351 if (!term) 4352 goto error; 4353 term->pow[poly->var] = i; 4354 term = isl_poly_foreach_term(rec->p[i], fn, term, user); 4355 if (!term) 4356 goto error; 4357 } 4358 term = isl_term_cow(term); 4359 if (!term) 4360 return NULL; 4361 term->pow[poly->var] = 0; 4362 4363 return term; 4364error: 4365 isl_term_free(term); 4366 return NULL; 4367} 4368 4369isl_stat isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp, 4370 isl_stat (*fn)(__isl_take isl_term *term, void *user), void *user) 4371{ 4372 isl_local *local; 4373 isl_term *term; 4374 4375 if (!qp) 4376 return isl_stat_error; 4377 4378 local = isl_qpolynomial_get_local(qp); 4379 term = isl_term_alloc(isl_space_copy(qp->dim), local); 4380 if (!term) 4381 return isl_stat_error; 4382 4383 term = isl_poly_foreach_term(qp->poly, fn, term, user); 4384 4385 isl_term_free(term); 4386 4387 return term ? isl_stat_ok : isl_stat_error; 4388} 4389 4390__isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term) 4391{ 4392 isl_poly *poly; 4393 isl_qpolynomial *qp; 4394 int i; 4395 isl_size n; 4396 4397 n = isl_term_dim(term, isl_dim_all); 4398 if (n < 0) 4399 term = isl_term_free(term); 4400 if (!term) 4401 return NULL; 4402 4403 poly = isl_poly_rat_cst(term->dim->ctx, term->n, term->d); 4404 for (i = 0; i < n; ++i) { 4405 if (!term->pow[i]) 4406 continue; 4407 poly = isl_poly_mul(poly, 4408 isl_poly_var_pow(term->dim->ctx, i, term->pow[i])); 4409 } 4410 4411 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), 4412 term->div->n_row, poly); 4413 if (!qp) 4414 goto error; 4415 isl_mat_free(qp->div); 4416 qp->div = isl_mat_copy(term->div); 4417 if (!qp->div) 4418 goto error; 4419 4420 isl_term_free(term); 4421 return qp; 4422error: 4423 isl_qpolynomial_free(qp); 4424 isl_term_free(term); 4425 return NULL; 4426} 4427 4428__isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp, 4429 __isl_take isl_space *space) 4430{ 4431 int i; 4432 int extra; 4433 isl_size total, d_set, d_qp; 4434 4435 if (!qp || !space) 4436 goto error; 4437 4438 if (isl_space_is_equal(qp->dim, space)) { 4439 isl_space_free(space); 4440 return qp; 4441 } 4442 4443 qp = isl_qpolynomial_cow(qp); 4444 if (!qp) 4445 goto error; 4446 4447 d_set = isl_space_dim(space, isl_dim_set); 4448 d_qp = isl_qpolynomial_domain_dim(qp, isl_dim_set); 4449 extra = d_set - d_qp; 4450 total = isl_space_dim(qp->dim, isl_dim_all); 4451 if (d_set < 0 || d_qp < 0 || total < 0) 4452 goto error; 4453 if (qp->div->n_row) { 4454 int *exp; 4455 4456 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row); 4457 if (!exp) 4458 goto error; 4459 for (i = 0; i < qp->div->n_row; ++i) 4460 exp[i] = extra + i; 4461 qp->poly = expand(qp->poly, exp, total); 4462 free(exp); 4463 if (!qp->poly) 4464 goto error; 4465 } 4466 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra); 4467 if (!qp->div) 4468 goto error; 4469 for (i = 0; i < qp->div->n_row; ++i) 4470 isl_seq_clr(qp->div->row[i] + 2 + total, extra); 4471 4472 isl_space_free(isl_qpolynomial_take_domain_space(qp)); 4473 qp = isl_qpolynomial_restore_domain_space(qp, space); 4474 4475 return qp; 4476error: 4477 isl_space_free(space); 4478 isl_qpolynomial_free(qp); 4479 return NULL; 4480} 4481 4482/* For each parameter or variable that does not appear in qp, 4483 * first eliminate the variable from all constraints and then set it to zero. 4484 */ 4485static __isl_give isl_set *fix_inactive(__isl_take isl_set *set, 4486 __isl_keep isl_qpolynomial *qp) 4487{ 4488 int *active = NULL; 4489 int i; 4490 isl_size d; 4491 isl_size nparam; 4492 isl_size nvar; 4493 4494 d = isl_set_dim(set, isl_dim_all); 4495 if (d < 0 || !qp) 4496 goto error; 4497 4498 active = isl_calloc_array(set->ctx, int, d); 4499 if (set_active(qp, active) < 0) 4500 goto error; 4501 4502 for (i = 0; i < d; ++i) 4503 if (!active[i]) 4504 break; 4505 4506 if (i == d) { 4507 free(active); 4508 return set; 4509 } 4510 4511 nparam = isl_set_dim(set, isl_dim_param); 4512 nvar = isl_set_dim(set, isl_dim_set); 4513 if (nparam < 0 || nvar < 0) 4514 goto error; 4515 for (i = 0; i < nparam; ++i) { 4516 if (active[i]) 4517 continue; 4518 set = isl_set_eliminate(set, isl_dim_param, i, 1); 4519 set = isl_set_fix_si(set, isl_dim_param, i, 0); 4520 } 4521 for (i = 0; i < nvar; ++i) { 4522 if (active[nparam + i]) 4523 continue; 4524 set = isl_set_eliminate(set, isl_dim_set, i, 1); 4525 set = isl_set_fix_si(set, isl_dim_set, i, 0); 4526 } 4527 4528 free(active); 4529 4530 return set; 4531error: 4532 free(active); 4533 isl_set_free(set); 4534 return NULL; 4535} 4536 4537struct isl_opt_data { 4538 isl_qpolynomial *qp; 4539 int first; 4540 isl_val *opt; 4541 int max; 4542}; 4543 4544static isl_stat opt_fn(__isl_take isl_point *pnt, void *user) 4545{ 4546 struct isl_opt_data *data = (struct isl_opt_data *)user; 4547 isl_val *val; 4548 4549 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt); 4550 if (data->first) { 4551 data->first = 0; 4552 data->opt = val; 4553 } else if (data->max) { 4554 data->opt = isl_val_max(data->opt, val); 4555 } else { 4556 data->opt = isl_val_min(data->opt, val); 4557 } 4558 4559 return isl_stat_ok; 4560} 4561 4562__isl_give isl_val *isl_qpolynomial_opt_on_domain( 4563 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max) 4564{ 4565 struct isl_opt_data data = { NULL, 1, NULL, max }; 4566 isl_bool is_cst; 4567 4568 if (!set || !qp) 4569 goto error; 4570 4571 is_cst = isl_poly_is_cst(qp->poly); 4572 if (is_cst < 0) 4573 goto error; 4574 if (is_cst) { 4575 isl_set_free(set); 4576 data.opt = isl_qpolynomial_get_constant_val(qp); 4577 isl_qpolynomial_free(qp); 4578 return data.opt; 4579 } 4580 4581 set = fix_inactive(set, qp); 4582 4583 data.qp = qp; 4584 if (isl_set_foreach_point(set, opt_fn, &data) < 0) 4585 goto error; 4586 4587 if (data.first) 4588 data.opt = isl_val_zero(isl_set_get_ctx(set)); 4589 4590 isl_set_free(set); 4591 isl_qpolynomial_free(qp); 4592 return data.opt; 4593error: 4594 isl_set_free(set); 4595 isl_qpolynomial_free(qp); 4596 isl_val_free(data.opt); 4597 return NULL; 4598} 4599 4600__isl_give isl_qpolynomial *isl_qpolynomial_morph_domain( 4601 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph) 4602{ 4603 int i; 4604 int n_sub; 4605 isl_ctx *ctx; 4606 isl_space *space; 4607 isl_poly **subs; 4608 isl_mat *mat, *diag; 4609 4610 qp = isl_qpolynomial_cow(qp); 4611 4612 space = isl_qpolynomial_peek_domain_space(qp); 4613 if (isl_morph_check_applies(morph, space) < 0) 4614 goto error; 4615 4616 ctx = isl_qpolynomial_get_ctx(qp); 4617 n_sub = morph->inv->n_row - 1; 4618 if (morph->inv->n_row != morph->inv->n_col) 4619 n_sub += qp->div->n_row; 4620 subs = isl_calloc_array(ctx, struct isl_poly *, n_sub); 4621 if (n_sub && !subs) 4622 goto error; 4623 4624 for (i = 0; 1 + i < morph->inv->n_row; ++i) 4625 subs[i] = isl_poly_from_affine(ctx, morph->inv->row[1 + i], 4626 morph->inv->row[0][0], morph->inv->n_col); 4627 if (morph->inv->n_row != morph->inv->n_col) 4628 for (i = 0; i < qp->div->n_row; ++i) 4629 subs[morph->inv->n_row - 1 + i] = 4630 isl_poly_var_pow(ctx, morph->inv->n_col - 1 + i, 1); 4631 4632 qp->poly = isl_poly_subs(qp->poly, 0, n_sub, subs); 4633 4634 for (i = 0; i < n_sub; ++i) 4635 isl_poly_free(subs[i]); 4636 free(subs); 4637 4638 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]); 4639 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv)); 4640 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]); 4641 mat = isl_mat_diagonal(mat, diag); 4642 qp->div = isl_mat_product(qp->div, mat); 4643 4644 if (!qp->poly || !qp->div) 4645 goto error; 4646 4647 isl_space_free(isl_qpolynomial_take_domain_space(qp)); 4648 space = isl_space_copy(morph->ran->dim); 4649 qp = isl_qpolynomial_restore_domain_space(qp, space); 4650 4651 isl_morph_free(morph); 4652 4653 return qp; 4654error: 4655 isl_qpolynomial_free(qp); 4656 isl_morph_free(morph); 4657 return NULL; 4658} 4659 4660__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul( 4661 __isl_take isl_union_pw_qpolynomial *upwqp1, 4662 __isl_take isl_union_pw_qpolynomial *upwqp2) 4663{ 4664 return isl_union_pw_qpolynomial_match_bin_op(upwqp1, upwqp2, 4665 &isl_pw_qpolynomial_mul); 4666} 4667 4668/* Reorder the dimension of "qp" according to the given reordering. 4669 */ 4670__isl_give isl_qpolynomial *isl_qpolynomial_realign_domain( 4671 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r) 4672{ 4673 isl_space *space; 4674 isl_poly *poly; 4675 isl_local *local; 4676 4677 if (!qp) 4678 goto error; 4679 4680 r = isl_reordering_extend(r, qp->div->n_row); 4681 if (!r) 4682 goto error; 4683 4684 local = isl_qpolynomial_take_local(qp); 4685 local = isl_local_reorder(local, isl_reordering_copy(r)); 4686 qp = isl_qpolynomial_restore_local(qp, local); 4687 4688 poly = isl_qpolynomial_take_poly(qp); 4689 poly = reorder(poly, r->pos); 4690 qp = isl_qpolynomial_restore_poly(qp, poly); 4691 4692 space = isl_reordering_get_space(r); 4693 qp = isl_qpolynomial_reset_domain_space(qp, space); 4694 4695 isl_reordering_free(r); 4696 return qp; 4697error: 4698 isl_qpolynomial_free(qp); 4699 isl_reordering_free(r); 4700 return NULL; 4701} 4702 4703__isl_give isl_qpolynomial *isl_qpolynomial_align_params( 4704 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model) 4705{ 4706 isl_space *domain_space; 4707 isl_bool equal_params; 4708 4709 domain_space = isl_qpolynomial_peek_domain_space(qp); 4710 equal_params = isl_space_has_equal_params(domain_space, model); 4711 if (equal_params < 0) 4712 goto error; 4713 if (!equal_params) { 4714 isl_reordering *exp; 4715 4716 exp = isl_parameter_alignment_reordering(domain_space, model); 4717 qp = isl_qpolynomial_realign_domain(qp, exp); 4718 } 4719 4720 isl_space_free(model); 4721 return qp; 4722error: 4723 isl_space_free(model); 4724 isl_qpolynomial_free(qp); 4725 return NULL; 4726} 4727 4728struct isl_split_periods_data { 4729 int max_periods; 4730 isl_pw_qpolynomial *res; 4731}; 4732 4733/* Create a slice where the integer division "div" has the fixed value "v". 4734 * In particular, if "div" refers to floor(f/m), then create a slice 4735 * 4736 * m v <= f <= m v + (m - 1) 4737 * 4738 * or 4739 * 4740 * f - m v >= 0 4741 * -f + m v + (m - 1) >= 0 4742 */ 4743static __isl_give isl_set *set_div_slice(__isl_take isl_space *space, 4744 __isl_keep isl_qpolynomial *qp, int div, isl_int v) 4745{ 4746 isl_size total; 4747 isl_basic_set *bset = NULL; 4748 int k; 4749 4750 total = isl_space_dim(space, isl_dim_all); 4751 if (total < 0 || !qp) 4752 goto error; 4753 4754 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0, 0, 2); 4755 4756 k = isl_basic_set_alloc_inequality(bset); 4757 if (k < 0) 4758 goto error; 4759 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total); 4760 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]); 4761 4762 k = isl_basic_set_alloc_inequality(bset); 4763 if (k < 0) 4764 goto error; 4765 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total); 4766 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]); 4767 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]); 4768 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1); 4769 4770 isl_space_free(space); 4771 return isl_set_from_basic_set(bset); 4772error: 4773 isl_basic_set_free(bset); 4774 isl_space_free(space); 4775 return NULL; 4776} 4777 4778static isl_stat split_periods(__isl_take isl_set *set, 4779 __isl_take isl_qpolynomial *qp, void *user); 4780 4781/* Create a slice of the domain "set" such that integer division "div" 4782 * has the fixed value "v" and add the results to data->res, 4783 * replacing the integer division by "v" in "qp". 4784 */ 4785static isl_stat set_div(__isl_take isl_set *set, 4786 __isl_take isl_qpolynomial *qp, int div, isl_int v, 4787 struct isl_split_periods_data *data) 4788{ 4789 int i; 4790 isl_size div_pos; 4791 isl_set *slice; 4792 isl_poly *cst; 4793 4794 slice = set_div_slice(isl_set_get_space(set), qp, div, v); 4795 set = isl_set_intersect(set, slice); 4796 4797 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 4798 if (div_pos < 0) 4799 goto error; 4800 4801 for (i = div + 1; i < qp->div->n_row; ++i) { 4802 if (isl_int_is_zero(qp->div->row[i][2 + div_pos + div])) 4803 continue; 4804 isl_int_addmul(qp->div->row[i][1], 4805 qp->div->row[i][2 + div_pos + div], v); 4806 isl_int_set_si(qp->div->row[i][2 + div_pos + div], 0); 4807 } 4808 4809 cst = isl_poly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one); 4810 qp = substitute_div(qp, div, cst); 4811 4812 return split_periods(set, qp, data); 4813error: 4814 isl_set_free(set); 4815 isl_qpolynomial_free(qp); 4816 return isl_stat_error; 4817} 4818 4819/* Split the domain "set" such that integer division "div" 4820 * has a fixed value (ranging from "min" to "max") on each slice 4821 * and add the results to data->res. 4822 */ 4823static isl_stat split_div(__isl_take isl_set *set, 4824 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max, 4825 struct isl_split_periods_data *data) 4826{ 4827 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) { 4828 isl_set *set_i = isl_set_copy(set); 4829 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp); 4830 4831 if (set_div(set_i, qp_i, div, min, data) < 0) 4832 goto error; 4833 } 4834 isl_set_free(set); 4835 isl_qpolynomial_free(qp); 4836 return isl_stat_ok; 4837error: 4838 isl_set_free(set); 4839 isl_qpolynomial_free(qp); 4840 return isl_stat_error; 4841} 4842 4843/* If "qp" refers to any integer division 4844 * that can only attain "max_periods" distinct values on "set" 4845 * then split the domain along those distinct values. 4846 * Add the results (or the original if no splitting occurs) 4847 * to data->res. 4848 */ 4849static isl_stat split_periods(__isl_take isl_set *set, 4850 __isl_take isl_qpolynomial *qp, void *user) 4851{ 4852 int i; 4853 isl_pw_qpolynomial *pwqp; 4854 struct isl_split_periods_data *data; 4855 isl_int min, max; 4856 isl_size div_pos; 4857 isl_stat r = isl_stat_ok; 4858 4859 data = (struct isl_split_periods_data *)user; 4860 4861 if (!set || !qp) 4862 goto error; 4863 4864 if (qp->div->n_row == 0) { 4865 pwqp = isl_pw_qpolynomial_alloc(set, qp); 4866 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp); 4867 return isl_stat_ok; 4868 } 4869 4870 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 4871 if (div_pos < 0) 4872 goto error; 4873 4874 isl_int_init(min); 4875 isl_int_init(max); 4876 for (i = 0; i < qp->div->n_row; ++i) { 4877 enum isl_lp_result lp_res; 4878 4879 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + div_pos, 4880 qp->div->n_row) != -1) 4881 continue; 4882 4883 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1, 4884 set->ctx->one, &min, NULL, NULL); 4885 if (lp_res == isl_lp_error) 4886 goto error2; 4887 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty) 4888 continue; 4889 isl_int_fdiv_q(min, min, qp->div->row[i][0]); 4890 4891 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1, 4892 set->ctx->one, &max, NULL, NULL); 4893 if (lp_res == isl_lp_error) 4894 goto error2; 4895 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty) 4896 continue; 4897 isl_int_fdiv_q(max, max, qp->div->row[i][0]); 4898 4899 isl_int_sub(max, max, min); 4900 if (isl_int_cmp_si(max, data->max_periods) < 0) { 4901 isl_int_add(max, max, min); 4902 break; 4903 } 4904 } 4905 4906 if (i < qp->div->n_row) { 4907 r = split_div(set, qp, i, min, max, data); 4908 } else { 4909 pwqp = isl_pw_qpolynomial_alloc(set, qp); 4910 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp); 4911 } 4912 4913 isl_int_clear(max); 4914 isl_int_clear(min); 4915 4916 return r; 4917error2: 4918 isl_int_clear(max); 4919 isl_int_clear(min); 4920error: 4921 isl_set_free(set); 4922 isl_qpolynomial_free(qp); 4923 return isl_stat_error; 4924} 4925 4926/* If any quasi-polynomial in pwqp refers to any integer division 4927 * that can only attain "max_periods" distinct values on its domain 4928 * then split the domain along those distinct values. 4929 */ 4930__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods( 4931 __isl_take isl_pw_qpolynomial *pwqp, int max_periods) 4932{ 4933 struct isl_split_periods_data data; 4934 4935 data.max_periods = max_periods; 4936 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp)); 4937 4938 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0) 4939 goto error; 4940 4941 isl_pw_qpolynomial_free(pwqp); 4942 4943 return data.res; 4944error: 4945 isl_pw_qpolynomial_free(data.res); 4946 isl_pw_qpolynomial_free(pwqp); 4947 return NULL; 4948} 4949 4950/* Construct a piecewise quasipolynomial that is constant on the given 4951 * domain. In particular, it is 4952 * 0 if cst == 0 4953 * 1 if cst == 1 4954 * infinity if cst == -1 4955 * 4956 * If cst == -1, then explicitly check whether the domain is empty and, 4957 * if so, return 0 instead. 4958 */ 4959static __isl_give isl_pw_qpolynomial *constant_on_domain( 4960 __isl_take isl_basic_set *bset, int cst) 4961{ 4962 isl_space *space; 4963 isl_qpolynomial *qp; 4964 4965 if (cst < 0 && isl_basic_set_is_empty(bset) == isl_bool_true) 4966 cst = 0; 4967 if (!bset) 4968 return NULL; 4969 4970 bset = isl_basic_set_params(bset); 4971 space = isl_basic_set_get_space(bset); 4972 if (cst < 0) 4973 qp = isl_qpolynomial_infty_on_domain(space); 4974 else if (cst == 0) 4975 qp = isl_qpolynomial_zero_on_domain(space); 4976 else 4977 qp = isl_qpolynomial_one_on_domain(space); 4978 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp); 4979} 4980 4981/* Internal data structure for multiplicative_call_factor_pw_qpolynomial. 4982 * "fn" is the function that is called on each factor. 4983 * "pwpq" collects the results. 4984 */ 4985struct isl_multiplicative_call_data_pw_qpolynomial { 4986 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset); 4987 isl_pw_qpolynomial *pwqp; 4988}; 4989 4990/* Call "fn" on "bset" and return the result, 4991 * but first check if "bset" has any redundant constraints or 4992 * implicit equality constraints. 4993 * If so, there may be further opportunities for detecting factors or 4994 * removing equality constraints, so recursively call 4995 * the top-level isl_basic_set_multiplicative_call. 4996 */ 4997static __isl_give isl_pw_qpolynomial *multiplicative_call_base( 4998 __isl_take isl_basic_set *bset, 4999 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset)) 5000{ 5001 isl_size n1, n2, n_eq; 5002 5003 n1 = isl_basic_set_n_constraint(bset); 5004 if (n1 < 0) 5005 bset = isl_basic_set_free(bset); 5006 bset = isl_basic_set_remove_redundancies(bset); 5007 bset = isl_basic_set_detect_equalities(bset); 5008 n2 = isl_basic_set_n_constraint(bset); 5009 n_eq = isl_basic_set_n_equality(bset); 5010 if (n2 < 0 || n_eq < 0) 5011 bset = isl_basic_set_free(bset); 5012 else if (n2 < n1 || n_eq > 0) 5013 return isl_basic_set_multiplicative_call(bset, fn); 5014 return fn(bset); 5015} 5016 5017/* isl_factorizer_every_factor_basic_set callback that applies 5018 * data->fn to the factor "bset" and multiplies in the result 5019 * in data->pwqp. 5020 */ 5021static isl_bool multiplicative_call_factor_pw_qpolynomial( 5022 __isl_keep isl_basic_set *bset, void *user) 5023{ 5024 struct isl_multiplicative_call_data_pw_qpolynomial *data = user; 5025 isl_pw_qpolynomial *res; 5026 5027 bset = isl_basic_set_copy(bset); 5028 res = multiplicative_call_base(bset, data->fn); 5029 data->pwqp = isl_pw_qpolynomial_mul(data->pwqp, res); 5030 if (!data->pwqp) 5031 return isl_bool_error; 5032 5033 return isl_bool_true; 5034} 5035 5036/* Factor bset, call fn on each of the factors and return the product. 5037 * 5038 * If no factors can be found, simply call fn on the input. 5039 * Otherwise, construct the factors based on the factorizer, 5040 * call fn on each factor and compute the product. 5041 */ 5042static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call( 5043 __isl_take isl_basic_set *bset, 5044 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset)) 5045{ 5046 struct isl_multiplicative_call_data_pw_qpolynomial data = { fn }; 5047 isl_space *space; 5048 isl_set *set; 5049 isl_factorizer *f; 5050 isl_qpolynomial *qp; 5051 isl_bool every; 5052 5053 f = isl_basic_set_factorizer(bset); 5054 if (!f) 5055 goto error; 5056 if (f->n_group == 0) { 5057 isl_factorizer_free(f); 5058 return multiplicative_call_base(bset, fn); 5059 } 5060 5061 space = isl_basic_set_get_space(bset); 5062 space = isl_space_params(space); 5063 set = isl_set_universe(isl_space_copy(space)); 5064 qp = isl_qpolynomial_one_on_domain(space); 5065 data.pwqp = isl_pw_qpolynomial_alloc(set, qp); 5066 5067 every = isl_factorizer_every_factor_basic_set(f, 5068 &multiplicative_call_factor_pw_qpolynomial, &data); 5069 if (every < 0) 5070 data.pwqp = isl_pw_qpolynomial_free(data.pwqp); 5071 5072 isl_basic_set_free(bset); 5073 isl_factorizer_free(f); 5074 5075 return data.pwqp; 5076error: 5077 isl_basic_set_free(bset); 5078 return NULL; 5079} 5080 5081/* Factor bset, call fn on each of the factors and return the product. 5082 * The function is assumed to evaluate to zero on empty domains, 5083 * to one on zero-dimensional domains and to infinity on unbounded domains 5084 * and will not be called explicitly on zero-dimensional or unbounded domains. 5085 * 5086 * We first check for some special cases and remove all equalities. 5087 * Then we hand over control to compressed_multiplicative_call. 5088 */ 5089__isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call( 5090 __isl_take isl_basic_set *bset, 5091 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset)) 5092{ 5093 isl_bool bounded; 5094 isl_size dim; 5095 isl_morph *morph; 5096 isl_pw_qpolynomial *pwqp; 5097 5098 if (!bset) 5099 return NULL; 5100 5101 if (isl_basic_set_plain_is_empty(bset)) 5102 return constant_on_domain(bset, 0); 5103 5104 dim = isl_basic_set_dim(bset, isl_dim_set); 5105 if (dim < 0) 5106 goto error; 5107 if (dim == 0) 5108 return constant_on_domain(bset, 1); 5109 5110 bounded = isl_basic_set_is_bounded(bset); 5111 if (bounded < 0) 5112 goto error; 5113 if (!bounded) 5114 return constant_on_domain(bset, -1); 5115 5116 if (bset->n_eq == 0) 5117 return compressed_multiplicative_call(bset, fn); 5118 5119 morph = isl_basic_set_full_compression(bset); 5120 bset = isl_morph_basic_set(isl_morph_copy(morph), bset); 5121 5122 pwqp = compressed_multiplicative_call(bset, fn); 5123 5124 morph = isl_morph_dom_params(morph); 5125 morph = isl_morph_ran_params(morph); 5126 morph = isl_morph_inverse(morph); 5127 5128 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph); 5129 5130 return pwqp; 5131error: 5132 isl_basic_set_free(bset); 5133 return NULL; 5134} 5135 5136/* Drop all floors in "qp", turning each integer division [a/m] into 5137 * a rational division a/m. If "down" is set, then the integer division 5138 * is replaced by (a-(m-1))/m instead. 5139 */ 5140static __isl_give isl_qpolynomial *qp_drop_floors( 5141 __isl_take isl_qpolynomial *qp, int down) 5142{ 5143 int i; 5144 isl_poly *s; 5145 5146 if (!qp) 5147 return NULL; 5148 if (qp->div->n_row == 0) 5149 return qp; 5150 5151 qp = isl_qpolynomial_cow(qp); 5152 if (!qp) 5153 return NULL; 5154 5155 for (i = qp->div->n_row - 1; i >= 0; --i) { 5156 if (down) { 5157 isl_int_sub(qp->div->row[i][1], 5158 qp->div->row[i][1], qp->div->row[i][0]); 5159 isl_int_add_ui(qp->div->row[i][1], 5160 qp->div->row[i][1], 1); 5161 } 5162 s = isl_poly_from_affine(qp->dim->ctx, qp->div->row[i] + 1, 5163 qp->div->row[i][0], qp->div->n_col - 1); 5164 qp = substitute_div(qp, i, s); 5165 if (!qp) 5166 return NULL; 5167 } 5168 5169 return qp; 5170} 5171 5172/* Drop all floors in "pwqp", turning each integer division [a/m] into 5173 * a rational division a/m. 5174 */ 5175static __isl_give isl_pw_qpolynomial *pwqp_drop_floors( 5176 __isl_take isl_pw_qpolynomial *pwqp) 5177{ 5178 int i; 5179 5180 if (!pwqp) 5181 return NULL; 5182 5183 if (isl_pw_qpolynomial_is_zero(pwqp)) 5184 return pwqp; 5185 5186 pwqp = isl_pw_qpolynomial_cow(pwqp); 5187 if (!pwqp) 5188 return NULL; 5189 5190 for (i = 0; i < pwqp->n; ++i) { 5191 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0); 5192 if (!pwqp->p[i].qp) 5193 goto error; 5194 } 5195 5196 return pwqp; 5197error: 5198 isl_pw_qpolynomial_free(pwqp); 5199 return NULL; 5200} 5201 5202/* Adjust all the integer divisions in "qp" such that they are at least 5203 * one over the given orthant (identified by "signs"). This ensures 5204 * that they will still be non-negative even after subtracting (m-1)/m. 5205 * 5206 * In particular, f is replaced by f' + v, changing f = [a/m] 5207 * to f' = [(a - m v)/m]. 5208 * If the constant term k in a is smaller than m, 5209 * the constant term of v is set to floor(k/m) - 1. 5210 * For any other term, if the coefficient c and the variable x have 5211 * the same sign, then no changes are needed. 5212 * Otherwise, if the variable is positive (and c is negative), 5213 * then the coefficient of x in v is set to floor(c/m). 5214 * If the variable is negative (and c is positive), 5215 * then the coefficient of x in v is set to ceil(c/m). 5216 */ 5217static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp, 5218 int *signs) 5219{ 5220 int i, j; 5221 isl_size div_pos; 5222 isl_vec *v = NULL; 5223 isl_poly *s; 5224 5225 qp = isl_qpolynomial_cow(qp); 5226 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div); 5227 if (div_pos < 0) 5228 return isl_qpolynomial_free(qp); 5229 qp->div = isl_mat_cow(qp->div); 5230 if (!qp->div) 5231 goto error; 5232 5233 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1); 5234 5235 for (i = 0; i < qp->div->n_row; ++i) { 5236 isl_int *row = qp->div->row[i]; 5237 v = isl_vec_clr(v); 5238 if (!v) 5239 goto error; 5240 if (isl_int_lt(row[1], row[0])) { 5241 isl_int_fdiv_q(v->el[0], row[1], row[0]); 5242 isl_int_sub_ui(v->el[0], v->el[0], 1); 5243 isl_int_submul(row[1], row[0], v->el[0]); 5244 } 5245 for (j = 0; j < div_pos; ++j) { 5246 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0) 5247 continue; 5248 if (signs[j] < 0) 5249 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]); 5250 else 5251 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]); 5252 isl_int_submul(row[2 + j], row[0], v->el[1 + j]); 5253 } 5254 for (j = 0; j < i; ++j) { 5255 if (isl_int_sgn(row[2 + div_pos + j]) >= 0) 5256 continue; 5257 isl_int_fdiv_q(v->el[1 + div_pos + j], 5258 row[2 + div_pos + j], row[0]); 5259 isl_int_submul(row[2 + div_pos + j], 5260 row[0], v->el[1 + div_pos + j]); 5261 } 5262 for (j = i + 1; j < qp->div->n_row; ++j) { 5263 if (isl_int_is_zero(qp->div->row[j][2 + div_pos + i])) 5264 continue; 5265 isl_seq_combine(qp->div->row[j] + 1, 5266 qp->div->ctx->one, qp->div->row[j] + 1, 5267 qp->div->row[j][2 + div_pos + i], v->el, 5268 v->size); 5269 } 5270 isl_int_set_si(v->el[1 + div_pos + i], 1); 5271 s = isl_poly_from_affine(qp->dim->ctx, v->el, 5272 qp->div->ctx->one, v->size); 5273 qp->poly = isl_poly_subs(qp->poly, div_pos + i, 1, &s); 5274 isl_poly_free(s); 5275 if (!qp->poly) 5276 goto error; 5277 } 5278 5279 isl_vec_free(v); 5280 return qp; 5281error: 5282 isl_vec_free(v); 5283 isl_qpolynomial_free(qp); 5284 return NULL; 5285} 5286 5287struct isl_to_poly_data { 5288 int sign; 5289 isl_pw_qpolynomial *res; 5290 isl_qpolynomial *qp; 5291}; 5292 5293/* Appoximate data->qp by a polynomial on the orthant identified by "signs". 5294 * We first make all integer divisions positive and then split the 5295 * quasipolynomials into terms with sign data->sign (the direction 5296 * of the requested approximation) and terms with the opposite sign. 5297 * In the first set of terms, each integer division [a/m] is 5298 * overapproximated by a/m, while in the second it is underapproximated 5299 * by (a-(m-1))/m. 5300 */ 5301static isl_stat to_polynomial_on_orthant(__isl_take isl_set *orthant, 5302 int *signs, void *user) 5303{ 5304 struct isl_to_poly_data *data = user; 5305 isl_pw_qpolynomial *t; 5306 isl_qpolynomial *qp, *up, *down; 5307 5308 qp = isl_qpolynomial_copy(data->qp); 5309 qp = make_divs_pos(qp, signs); 5310 5311 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign); 5312 up = qp_drop_floors(up, 0); 5313 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign); 5314 down = qp_drop_floors(down, 1); 5315 5316 isl_qpolynomial_free(qp); 5317 qp = isl_qpolynomial_add(up, down); 5318 5319 t = isl_pw_qpolynomial_alloc(orthant, qp); 5320 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t); 5321 5322 return isl_stat_ok; 5323} 5324 5325/* Approximate each quasipolynomial by a polynomial. If "sign" is positive, 5326 * the polynomial will be an overapproximation. If "sign" is negative, 5327 * it will be an underapproximation. If "sign" is zero, the approximation 5328 * will lie somewhere in between. 5329 * 5330 * In particular, is sign == 0, we simply drop the floors, turning 5331 * the integer divisions into rational divisions. 5332 * Otherwise, we split the domains into orthants, make all integer divisions 5333 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m, 5334 * depending on the requested sign and the sign of the term in which 5335 * the integer division appears. 5336 */ 5337__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial( 5338 __isl_take isl_pw_qpolynomial *pwqp, int sign) 5339{ 5340 int i; 5341 struct isl_to_poly_data data; 5342 5343 if (sign == 0) 5344 return pwqp_drop_floors(pwqp); 5345 5346 if (!pwqp) 5347 return NULL; 5348 5349 data.sign = sign; 5350 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp)); 5351 5352 for (i = 0; i < pwqp->n; ++i) { 5353 if (pwqp->p[i].qp->div->n_row == 0) { 5354 isl_pw_qpolynomial *t; 5355 t = isl_pw_qpolynomial_alloc( 5356 isl_set_copy(pwqp->p[i].set), 5357 isl_qpolynomial_copy(pwqp->p[i].qp)); 5358 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t); 5359 continue; 5360 } 5361 data.qp = pwqp->p[i].qp; 5362 if (isl_set_foreach_orthant(pwqp->p[i].set, 5363 &to_polynomial_on_orthant, &data) < 0) 5364 goto error; 5365 } 5366 5367 isl_pw_qpolynomial_free(pwqp); 5368 5369 return data.res; 5370error: 5371 isl_pw_qpolynomial_free(pwqp); 5372 isl_pw_qpolynomial_free(data.res); 5373 return NULL; 5374} 5375 5376static __isl_give isl_pw_qpolynomial *poly_entry( 5377 __isl_take isl_pw_qpolynomial *pwqp, void *user) 5378{ 5379 int *sign = user; 5380 5381 return isl_pw_qpolynomial_to_polynomial(pwqp, *sign); 5382} 5383 5384__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial( 5385 __isl_take isl_union_pw_qpolynomial *upwqp, int sign) 5386{ 5387 return isl_union_pw_qpolynomial_transform_inplace(upwqp, 5388 &poly_entry, &sign); 5389} 5390 5391__isl_give isl_basic_map *isl_basic_map_from_qpolynomial( 5392 __isl_take isl_qpolynomial *qp) 5393{ 5394 isl_local_space *ls; 5395 isl_vec *vec; 5396 isl_aff *aff; 5397 isl_basic_map *bmap; 5398 isl_bool is_affine; 5399 5400 if (!qp) 5401 return NULL; 5402 is_affine = isl_poly_is_affine(qp->poly); 5403 if (is_affine < 0) 5404 goto error; 5405 if (!is_affine) 5406 isl_die(qp->dim->ctx, isl_error_invalid, 5407 "input quasi-polynomial not affine", goto error); 5408 ls = isl_qpolynomial_get_domain_local_space(qp); 5409 vec = isl_qpolynomial_extract_affine(qp); 5410 aff = isl_aff_alloc_vec(ls, vec); 5411 bmap = isl_basic_map_from_aff(aff); 5412 isl_qpolynomial_free(qp); 5413 return bmap; 5414error: 5415 isl_qpolynomial_free(qp); 5416 return NULL; 5417} 5418