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