isl_polynomial.c

   1 #include <isl_seq.h>
   2 #include <isl_polynomial_private.h>
   3 #include <isl_point_private.h>
   4 #include <isl_dim_private.h>
   5 #include <isl_map_private.h>
   6
   7 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
   8 {
   9         if (!up)
  10                 return -1;
  11
  12         return up->var < 0;
  13 }
  14
  15 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
  16 {
  17         if (!up)
  18                 return NULL;
  19
  20         isl_assert(up->ctx, up->var < 0, return NULL);
  21
  22         return (struct isl_upoly_cst *)up;
  23 }
  24
  25 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
  26 {
  27         if (!up)
  28                 return NULL;
  29
  30         isl_assert(up->ctx, up->var >= 0, return NULL);
  31
  32         return (struct isl_upoly_rec *)up;
  33 }
  34
  35 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
  36         __isl_keep struct isl_upoly *up2)
  37 {
  38         int i;
  39         struct isl_upoly_rec *rec1, *rec2;
  40
  41         if (!up1 || !up2)
  42                 return -1;
  43         if (up1 == up2)
  44                 return 1;
  45         if (up1->var != up2->var)
  46                 return 0;
  47         if (isl_upoly_is_cst(up1)) {
  48                 struct isl_upoly_cst *cst1, *cst2;
  49                 cst1 = isl_upoly_as_cst(up1);
  50                 cst2 = isl_upoly_as_cst(up2);
  51                 if (!cst1 || !cst2)
  52                         return -1;
  53                 return isl_int_eq(cst1->n, cst2->n) &&
  54                        isl_int_eq(cst1->d, cst2->d);
  55         }
  56
  57         rec1 = isl_upoly_as_rec(up1);
  58         rec2 = isl_upoly_as_rec(up2);
  59         if (!rec1 || !rec2)
  60                 return -1;
  61
  62         if (rec1->n != rec2->n)
  63                 return 0;
  64
  65         for (i = 0; i < rec1->n; ++i) {
  66                 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
  67                 if (eq < 0 || !eq)
  68                         return eq;
  69         }
  70
  71         return 1;
  72 }
  73
  74 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
  75 {
  76         struct isl_upoly_cst *cst;
  77
  78         if (!up)
  79                 return -1;
  80         if (!isl_upoly_is_cst(up))
  81                 return 0;
  82
  83         cst = isl_upoly_as_cst(up);
  84         if (!cst)
  85                 return -1;
  86
  87         return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
  88 }
  89
  90 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
  91 {
  92         struct isl_upoly_cst *cst;
  93
  94         if (!up)
  95                 return -1;
  96         if (!isl_upoly_is_cst(up))
  97                 return 0;
  98
  99         cst = isl_upoly_as_cst(up);
 100         if (!cst)
 101                 return -1;
 102
 103         return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
 104 }
 105
 106 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
 107 {
 108         struct isl_upoly_cst *cst;
 109
 110         if (!up)
 111                 return -1;
 112         if (!isl_upoly_is_cst(up))
 113                 return 0;
 114
 115         cst = isl_upoly_as_cst(up);
 116         if (!cst)
 117                 return -1;
 118
 119         return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
 120 }
 121
 122 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
 123 {
 124         struct isl_upoly_cst *cst;
 125
 126         if (!up)
 127                 return -1;
 128         if (!isl_upoly_is_cst(up))
 129                 return 0;
 130
 131         cst = isl_upoly_as_cst(up);
 132         if (!cst)
 133                 return -1;
 134
 135         return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
 136 }
 137
 138 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
 139 {
 140         struct isl_upoly_cst *cst;
 141
 142         if (!up)
 143                 return -1;
 144         if (!isl_upoly_is_cst(up))
 145                 return 0;
 146
 147         cst = isl_upoly_as_cst(up);
 148         if (!cst)
 149                 return -1;
 150
 151         return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
 152 }
 153
 154 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
 155 {
 156         struct isl_upoly_cst *cst;
 157
 158         if (!up)
 159                 return -1;
 160         if (!isl_upoly_is_cst(up))
 161                 return 0;
 162
 163         cst = isl_upoly_as_cst(up);
 164         if (!cst)
 165                 return -1;
 166
 167         return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
 168 }
 169
 170 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
 171 {
 172         struct isl_upoly_cst *cst;
 173
 174         cst = isl_alloc_type(ctx, struct isl_upoly_cst);
 175         if (!cst)
 176                 return NULL;
 177
 178         cst->up.ref = 1;
 179         cst->up.ctx = ctx;
 180         isl_ctx_ref(ctx);
 181         cst->up.var = -1;
 182
 183         isl_int_init(cst->n);
 184         isl_int_init(cst->d);
 185
 186         return cst;
 187 }
 188
 189 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
 190 {
 191         struct isl_upoly_cst *cst;
 192
 193         cst = isl_upoly_cst_alloc(ctx);
 194         if (!cst)
 195                 return NULL;
 196
 197         isl_int_set_si(cst->n, 0);
 198         isl_int_set_si(cst->d, 1);
 199
 200         return &cst->up;
 201 }
 202
 203 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
 204 {
 205         struct isl_upoly_cst *cst;
 206
 207         cst = isl_upoly_cst_alloc(ctx);
 208         if (!cst)
 209                 return NULL;
 210
 211         isl_int_set_si(cst->n, 1);
 212         isl_int_set_si(cst->d, 0);
 213
 214         return &cst->up;
 215 }
 216
 217 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
 218 {
 219         struct isl_upoly_cst *cst;
 220
 221         cst = isl_upoly_cst_alloc(ctx);
 222         if (!cst)
 223                 return NULL;
 224
 225         isl_int_set_si(cst->n, 0);
 226         isl_int_set_si(cst->d, 0);
 227
 228         return &cst->up;
 229 }
 230
 231 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
 232         isl_int n, isl_int d)
 233 {
 234         struct isl_upoly_cst *cst;
 235
 236         cst = isl_upoly_cst_alloc(ctx);
 237         if (!cst)
 238                 return NULL;
 239
 240         isl_int_set(cst->n, n);
 241         isl_int_set(cst->d, d);
 242
 243         return &cst->up;
 244 }
 245
 246 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
 247         int var, int size)
 248 {
 249         struct isl_upoly_rec *rec;
 250
 251         isl_assert(ctx, var >= 0, return NULL);
 252         isl_assert(ctx, size >= 0, return NULL);
 253         rec = isl_calloc(dim->ctx, struct isl_upoly_rec,
 254                         sizeof(struct isl_upoly_rec) +
 255                         (size - 1) * sizeof(struct isl_upoly *));
 256         if (!rec)
 257                 return NULL;
 258
 259         rec->up.ref = 1;
 260         rec->up.ctx = ctx;
 261         isl_ctx_ref(ctx);
 262         rec->up.var = var;
 263
 264         rec->n = 0;
 265         rec->size = size;
 266
 267         return rec;
 268 error:
 269         isl_upoly_free(&rec->up);
 270         return NULL;
 271 }
 272
 273 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
 274 {
 275         struct isl_upoly_cst *cst;
 276
 277         if (!qp)
 278                 return -1;
 279
 280         return isl_upoly_is_zero(qp->upoly);
 281 }
 282
 283 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
 284 {
 285         struct isl_upoly_cst *cst;
 286
 287         if (!qp)
 288                 return -1;
 289
 290         return isl_upoly_is_one(qp->upoly);
 291 }
 292
 293 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
 294 {
 295         isl_int_clear(cst->n);
 296         isl_int_clear(cst->d);
 297 }
 298
 299 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
 300 {
 301         int i;
 302
 303         for (i = 0; i < rec->n; ++i)
 304                 isl_upoly_free(rec->p[i]);
 305 }
 306
 307 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
 308 {
 309         if (!up)
 310                 return NULL;
 311
 312         up->ref++;
 313         return up;
 314 }
 315
 316 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
 317 {
 318         struct isl_upoly_cst *cst;
 319         struct isl_upoly_cst *dup;
 320
 321         cst = isl_upoly_as_cst(up);
 322         if (!cst)
 323                 return NULL;
 324
 325         dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
 326         if (!dup)
 327                 return NULL;
 328         isl_int_set(dup->n, cst->n);
 329         isl_int_set(dup->d, cst->d);
 330
 331         return &dup->up;
 332 }
 333
 334 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
 335 {
 336         int i;
 337         struct isl_upoly_rec *rec;
 338         struct isl_upoly_rec *dup;
 339
 340         rec = isl_upoly_as_rec(up);
 341         if (!rec)
 342                 return NULL;
 343
 344         dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
 345         if (!dup)
 346                 return NULL;
 347
 348         for (i = 0; i < rec->n; ++i) {
 349                 dup->p[i] = isl_upoly_copy(rec->p[i]);
 350                 if (!dup->p[i])
 351                         goto error;
 352                 dup->n++;
 353         }
 354
 355         return &dup->up;
 356 error:
 357         isl_upoly_free(&dup->up);
 358         return NULL;
 359 }
 360
 361 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
 362 {
 363         struct isl_upoly *dup;
 364
 365         if (!up)
 366                 return NULL;
 367
 368         if (isl_upoly_is_cst(up))
 369                 return isl_upoly_dup_cst(up);
 370         else
 371                 return isl_upoly_dup_rec(up);
 372 }
 373
 374 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
 375 {
 376         if (!up)
 377                 return NULL;
 378
 379         if (up->ref == 1)
 380                 return up;
 381         up->ref--;
 382         return isl_upoly_dup(up);
 383 }
 384
 385 void isl_upoly_free(__isl_take struct isl_upoly *up)
 386 {
 387         if (!up)
 388                 return;
 389
 390         if (--up->ref > 0)
 391                 return;
 392
 393         if (up->var < 0)
 394                 upoly_free_cst((struct isl_upoly_cst *)up);
 395         else
 396                 upoly_free_rec((struct isl_upoly_rec *)up);
 397
 398         isl_ctx_deref(up->ctx);
 399         free(up);
 400 }
 401
 402 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
 403 {
 404         isl_int gcd;
 405
 406         isl_int_init(gcd);
 407         isl_int_gcd(gcd, cst->n, cst->d);
 408         if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
 409                 isl_int_divexact(cst->n, cst->n, gcd);
 410                 isl_int_divexact(cst->d, cst->d, gcd);
 411         }
 412         isl_int_clear(gcd);
 413 }
 414
 415 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
 416         __isl_take struct isl_upoly *up2)
 417 {
 418         struct isl_upoly_cst *cst1;
 419         struct isl_upoly_cst *cst2;
 420
 421         up1 = isl_upoly_cow(up1);
 422         if (!up1 || !up2)
 423                 goto error;
 424
 425         cst1 = isl_upoly_as_cst(up1);
 426         cst2 = isl_upoly_as_cst(up2);
 427
 428         if (isl_int_eq(cst1->d, cst2->d))
 429                 isl_int_add(cst1->n, cst1->n, cst2->n);
 430         else {
 431                 isl_int_mul(cst1->n, cst1->n, cst2->d);
 432                 isl_int_addmul(cst1->n, cst2->n, cst1->d);
 433                 isl_int_mul(cst1->d, cst1->d, cst2->d);
 434         }
 435
 436         isl_upoly_cst_reduce(cst1);
 437
 438         isl_upoly_free(up2);
 439         return up1;
 440 error:
 441         isl_upoly_free(up1);
 442         isl_upoly_free(up2);
 443         return NULL;
 444 }
 445
 446 static __isl_give struct isl_upoly *replace_by_zero(
 447         __isl_take struct isl_upoly *up)
 448 {
 449         struct isl_ctx *ctx;
 450
 451         if (!up)
 452                 return NULL;
 453         ctx = up->ctx;
 454         isl_upoly_free(up);
 455         return isl_upoly_zero(ctx);
 456 }
 457
 458 static __isl_give struct isl_upoly *replace_by_constant_term(
 459         __isl_take struct isl_upoly *up)
 460 {
 461         struct isl_upoly_rec *rec;
 462         struct isl_upoly *cst;
 463
 464         if (!up)
 465                 return NULL;
 466
 467         rec = isl_upoly_as_rec(up);
 468         if (!rec)
 469                 goto error;
 470         isl_assert(up->ctx, rec->n == 1, goto error);
 471         cst = isl_upoly_copy(rec->p[0]);
 472         isl_upoly_free(up);
 473         return cst;
 474 error:
 475         isl_upoly_free(up);
 476         return NULL;
 477 }
 478
 479 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
 480         __isl_take struct isl_upoly *up2)
 481 {
 482         int i;
 483         struct isl_upoly_rec *rec1, *rec2;
 484
 485         if (!up1 || !up2)
 486                 goto error;
 487
 488         if (isl_upoly_is_nan(up1)) {
 489                 isl_upoly_free(up2);
 490                 return up1;
 491         }
 492
 493         if (isl_upoly_is_nan(up2)) {
 494                 isl_upoly_free(up1);
 495                 return up2;
 496         }
 497
 498         if (isl_upoly_is_zero(up1)) {
 499                 isl_upoly_free(up1);
 500                 return up2;
 501         }
 502
 503         if (isl_upoly_is_zero(up2)) {
 504                 isl_upoly_free(up2);
 505                 return up1;
 506         }
 507
 508         if (up1->var < up2->var)
 509                 return isl_upoly_sum(up2, up1);
 510
 511         if (up2->var < up1->var) {
 512                 struct isl_upoly_rec *rec;
 513                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 514                         isl_upoly_free(up1);
 515                         return up2;
 516                 }
 517                 up1 = isl_upoly_cow(up1);
 518                 rec = isl_upoly_as_rec(up1);
 519                 if (!rec)
 520                         goto error;
 521                 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
 522                 if (rec->n == 1)
 523                         up1 = replace_by_constant_term(up1);
 524                 return up1;
 525         }
 526
 527         if (isl_upoly_is_cst(up1))
 528                 return isl_upoly_sum_cst(up1, up2);
 529
 530         rec1 = isl_upoly_as_rec(up1);
 531         rec2 = isl_upoly_as_rec(up2);
 532         if (!rec1 || !rec2)
 533                 goto error;
 534
 535         if (rec1->n < rec2->n)
 536                 return isl_upoly_sum(up2, up1);
 537
 538         up1 = isl_upoly_cow(up1);
 539         rec1 = isl_upoly_as_rec(up1);
 540         if (!rec1)
 541                 goto error;
 542
 543         for (i = rec2->n - 1; i >= 0; --i) {
 544                 rec1->p[i] = isl_upoly_sum(rec1->p[i],
 545                                             isl_upoly_copy(rec2->p[i]));
 546                 if (!rec1->p[i])
 547                         goto error;
 548                 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
 549                         isl_upoly_free(rec1->p[i]);
 550                         rec1->n--;
 551                 }
 552         }
 553
 554         if (rec1->n == 0)
 555                 up1 = replace_by_zero(up1);
 556         else if (rec1->n == 1)
 557                 up1 = replace_by_constant_term(up1);
 558
 559         isl_upoly_free(up2);
 560
 561         return up1;
 562 error:
 563         isl_upoly_free(up1);
 564         isl_upoly_free(up2);
 565         return NULL;
 566 }
 567
 568 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
 569 {
 570         struct isl_upoly_cst *cst;
 571
 572         if (isl_upoly_is_zero(up))
 573                 return up;
 574
 575         up = isl_upoly_cow(up);
 576         if (!up)
 577                 return NULL;
 578
 579         cst = isl_upoly_as_cst(up);
 580
 581         isl_int_neg(cst->n, cst->n);
 582
 583         return up;
 584 }
 585
 586 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
 587 {
 588         int i;
 589         struct isl_upoly_rec *rec;
 590
 591         if (!up)
 592                 return NULL;
 593
 594         if (isl_upoly_is_cst(up))
 595                 return isl_upoly_neg_cst(up);
 596
 597         up = isl_upoly_cow(up);
 598         rec = isl_upoly_as_rec(up);
 599         if (!rec)
 600                 goto error;
 601
 602         for (i = 0; i < rec->n; ++i) {
 603                 rec->p[i] = isl_upoly_neg(rec->p[i]);
 604                 if (!rec->p[i])
 605                         goto error;
 606         }
 607
 608         return up;
 609 error:
 610         isl_upoly_free(up);
 611         return NULL;
 612 }
 613
 614 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
 615         __isl_take struct isl_upoly *up2)
 616 {
 617         struct isl_upoly_cst *cst1;
 618         struct isl_upoly_cst *cst2;
 619
 620         up1 = isl_upoly_cow(up1);
 621         if (!up1 || !up2)
 622                 goto error;
 623
 624         cst1 = isl_upoly_as_cst(up1);
 625         cst2 = isl_upoly_as_cst(up2);
 626
 627         isl_int_mul(cst1->n, cst1->n, cst2->n);
 628         isl_int_mul(cst1->d, cst1->d, cst2->d);
 629
 630         isl_upoly_cst_reduce(cst1);
 631
 632         isl_upoly_free(up2);
 633         return up1;
 634 error:
 635         isl_upoly_free(up1);
 636         isl_upoly_free(up2);
 637         return NULL;
 638 }
 639
 640 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
 641         __isl_take struct isl_upoly *up2)
 642 {
 643         struct isl_upoly_rec *rec1;
 644         struct isl_upoly_rec *rec2;
 645         struct isl_upoly_rec *res;
 646         int i, j;
 647         int size;
 648
 649         rec1 = isl_upoly_as_rec(up1);
 650         rec2 = isl_upoly_as_rec(up2);
 651         if (!rec1 || !rec2)
 652                 goto error;
 653         size = rec1->n + rec2->n - 1;
 654         res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
 655         if (!res)
 656                 goto error;
 657
 658         for (i = 0; i < rec1->n; ++i) {
 659                 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
 660                                             isl_upoly_copy(rec1->p[i]));
 661                 if (!res->p[i])
 662                         goto error;
 663                 res->n++;
 664         }
 665         for (; i < size; ++i) {
 666                 res->p[i] = isl_upoly_zero(up1->ctx);
 667                 if (!res->p[i])
 668                         goto error;
 669                 res->n++;
 670         }
 671         for (i = 0; i < rec1->n; ++i) {
 672                 for (j = 1; j < rec2->n; ++j) {
 673                         struct isl_upoly *up;
 674                         up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
 675                                             isl_upoly_copy(rec1->p[i]));
 676                         res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
 677                         if (!res->p[i + j])
 678                                 goto error;
 679                 }
 680         }
 681
 682         isl_upoly_free(up1);
 683         isl_upoly_free(up2);
 684
 685         return &res->up;
 686 error:
 687         isl_upoly_free(up1);
 688         isl_upoly_free(up2);
 689         isl_upoly_free(&res->up);
 690         return NULL;
 691 }
 692
 693 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
 694         __isl_take struct isl_upoly *up2)
 695 {
 696         if (!up1 || !up2)
 697                 goto error;
 698
 699         if (isl_upoly_is_nan(up1)) {
 700                 isl_upoly_free(up2);
 701                 return up1;
 702         }
 703
 704         if (isl_upoly_is_nan(up2)) {
 705                 isl_upoly_free(up1);
 706                 return up2;
 707         }
 708
 709         if (isl_upoly_is_zero(up1)) {
 710                 isl_upoly_free(up2);
 711                 return up1;
 712         }
 713
 714         if (isl_upoly_is_zero(up2)) {
 715                 isl_upoly_free(up1);
 716                 return up2;
 717         }
 718
 719         if (isl_upoly_is_one(up1)) {
 720                 isl_upoly_free(up1);
 721                 return up2;
 722         }
 723
 724         if (isl_upoly_is_one(up2)) {
 725                 isl_upoly_free(up2);
 726                 return up1;
 727         }
 728
 729         if (up1->var < up2->var)
 730                 return isl_upoly_mul(up2, up1);
 731
 732         if (up2->var < up1->var) {
 733                 int i;
 734                 struct isl_upoly_rec *rec;
 735                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 736                         isl_ctx *ctx = up1->ctx;
 737                         isl_upoly_free(up1);
 738                         isl_upoly_free(up2);
 739                         return isl_upoly_nan(ctx);
 740                 }
 741                 up1 = isl_upoly_cow(up1);
 742                 rec = isl_upoly_as_rec(up1);
 743                 if (!rec)
 744                         goto error;
 745
 746                 for (i = 0; i < rec->n; ++i) {
 747                         rec->p[i] = isl_upoly_mul(rec->p[i],
 748                                                     isl_upoly_copy(up2));
 749                         if (!rec->p[i])
 750                                 goto error;
 751                 }
 752                 isl_upoly_free(up2);
 753                 return up1;
 754         }
 755
 756         if (isl_upoly_is_cst(up1))
 757                 return isl_upoly_mul_cst(up1, up2);
 758
 759         return isl_upoly_mul_rec(up1, up2);
 760 error:
 761         isl_upoly_free(up1);
 762         isl_upoly_free(up2);
 763         return NULL;
 764 }
 765
 766 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
 767         unsigned n_div)
 768 {
 769         struct isl_qpolynomial *qp;
 770         unsigned total;
 771
 772         if (!dim)
 773                 return NULL;
 774
 775         total = isl_dim_total(dim);
 776
 777         qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
 778         if (!qp)
 779                 return NULL;
 780
 781         qp->ref = 1;
 782         qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
 783         if (!qp->div)
 784                 goto error;
 785
 786         qp->dim = dim;
 787
 788         return qp;
 789 error:
 790         isl_dim_free(dim);
 791         isl_qpolynomial_free(qp);
 792         return NULL;
 793 }
 794
 795 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
 796 {
 797         if (!qp)
 798                 return NULL;
 799
 800         qp->ref++;
 801         return qp;
 802 }
 803
 804 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
 805 {
 806         struct isl_qpolynomial *dup;
 807
 808         if (!qp)
 809                 return NULL;
 810
 811         dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row);
 812         if (!dup)
 813                 return NULL;
 814         isl_mat_free(dup->div);
 815         dup->div = isl_mat_copy(qp->div);
 816         if (!dup->div)
 817                 goto error;
 818         dup->upoly = isl_upoly_copy(qp->upoly);
 819         if (!dup->upoly)
 820                 goto error;
 821
 822         return dup;
 823 error:
 824         isl_qpolynomial_free(dup);
 825         return NULL;
 826 }
 827
 828 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
 829 {
 830         if (!qp)
 831                 return NULL;
 832
 833         if (qp->ref == 1)
 834                 return qp;
 835         qp->ref--;
 836         return isl_qpolynomial_dup(qp);
 837 }
 838
 839 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
 840 {
 841         if (!qp)
 842                 return;
 843
 844         if (--qp->ref > 0)
 845                 return;
 846
 847         isl_dim_free(qp->dim);
 848         isl_mat_free(qp->div);
 849         isl_upoly_free(qp->upoly);
 850
 851         free(qp);
 852 }
 853
 854 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
 855 {
 856         int n_row, n_col;
 857         int equal;
 858
 859         isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
 860                                 div1->n_col >= div2->n_col, return -1);
 861
 862         if (div1->n_row == div2->n_row)
 863                 return isl_mat_is_equal(div1, div2);
 864
 865         n_row = div1->n_row;
 866         n_col = div1->n_col;
 867         div1->n_row = div2->n_row;
 868         div1->n_col = div2->n_col;
 869
 870         equal = isl_mat_is_equal(div1, div2);
 871
 872         div1->n_row = n_row;
 873         div1->n_col = n_col;
 874
 875         return equal;
 876 }
 877
 878 static void expand_row(__isl_keep isl_mat *dst, int d,
 879         __isl_keep isl_mat *src, int s, int *exp)
 880 {
 881         int i;
 882         unsigned c = src->n_col - src->n_row;
 883
 884         isl_seq_cpy(dst->row[d], src->row[s], c);
 885         isl_seq_clr(dst->row[d] + c, dst->n_col - c);
 886
 887         for (i = 0; i < s; ++i)
 888                 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
 889 }
 890
 891 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
 892 {
 893         int li, lj;
 894
 895         li = isl_seq_last_non_zero(div->row[i], div->n_col);
 896         lj = isl_seq_last_non_zero(div->row[j], div->n_col);
 897
 898         if (li != lj)
 899                 return li - lj;
 900
 901         return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
 902 }
 903
 904 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
 905         __isl_keep isl_mat *div2, int *exp1, int *exp2)
 906 {
 907         int i, j, k;
 908         isl_mat *div = NULL;
 909         unsigned d = div1->n_col - div1->n_row;
 910
 911         div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
 912                                 d + div1->n_row + div2->n_row);
 913         if (!div)
 914                 return NULL;
 915
 916         for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
 917                 int cmp;
 918
 919                 expand_row(div, k, div1, i, exp1);
 920                 expand_row(div, k + 1, div2, j, exp2);
 921
 922                 cmp = cmp_row(div, k, k + 1);
 923                 if (cmp == 0) {
 924                         exp1[i++] = k;
 925                         exp2[j++] = k;
 926                 } else if (cmp < 0) {
 927                         exp1[i++] = k;
 928                 } else {
 929                         exp2[j++] = k;
 930                         isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
 931                 }
 932         }
 933         for (; i < div1->n_row; ++i, ++k) {
 934                 expand_row(div, k, div1, i, exp1);
 935                 exp1[i] = k;
 936         }
 937         for (; j < div2->n_row; ++j, ++k) {
 938                 expand_row(div, k, div2, j, exp2);
 939                 exp2[j] = k;
 940         }
 941
 942         div->n_row = k;
 943         div->n_col = d + k;
 944
 945         return div;
 946 }
 947
 948 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
 949         int *exp, int first)
 950 {
 951         int i;
 952         struct isl_upoly_rec *rec;
 953
 954         if (isl_upoly_is_cst(up))
 955                 return up;
 956
 957         if (up->var < first)
 958                 return up;
 959
 960         if (exp[up->var - first] == up->var - first)
 961                 return up;
 962
 963         up = isl_upoly_cow(up);
 964         if (!up)
 965                 goto error;
 966
 967         up->var = exp[up->var - first] + first;
 968
 969         rec = isl_upoly_as_rec(up);
 970         if (!rec)
 971                 goto error;
 972
 973         for (i = 0; i < rec->n; ++i) {
 974                 rec->p[i] = expand(rec->p[i], exp, first);
 975                 if (!rec->p[i])
 976                         goto error;
 977         }
 978
 979         return up;
 980 error:
 981         isl_upoly_free(up);
 982         return NULL;
 983 }
 984
 985 static __isl_give isl_qpolynomial *with_merged_divs(
 986         __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
 987                                           __isl_take isl_qpolynomial *qp2),
 988         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
 989 {
 990         int *exp1 = NULL;
 991         int *exp2 = NULL;
 992         isl_mat *div = NULL;
 993
 994         qp1 = isl_qpolynomial_cow(qp1);
 995         qp2 = isl_qpolynomial_cow(qp2);
 996
 997         if (!qp1 || !qp2)
 998                 goto error;
 999
1000         isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1001                                 qp1->div->n_col >= qp2->div->n_col, goto error);
1002
1003         exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1004         exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1005         if (!exp1 || !exp2)
1006                 goto error;
1007
1008         div = merge_divs(qp1->div, qp2->div, exp1, exp2);
1009         if (!div)
1010                 goto error;
1011
1012         isl_mat_free(qp1->div);
1013         qp1->div = isl_mat_copy(div);
1014         isl_mat_free(qp2->div);
1015         qp2->div = isl_mat_copy(div);
1016
1017         qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1018         qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1019
1020         if (!qp1->upoly || !qp2->upoly)
1021                 goto error;
1022
1023         isl_mat_free(div);
1024         free(exp1);
1025         free(exp2);
1026
1027         return fn(qp1, qp2);
1028 error:
1029         isl_mat_free(div);
1030         free(exp1);
1031         free(exp2);
1032         isl_qpolynomial_free(qp1);
1033         isl_qpolynomial_free(qp2);
1034         return NULL;
1035 }
1036
1037 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1038         __isl_take isl_qpolynomial *qp2)
1039 {
1040         qp1 = isl_qpolynomial_cow(qp1);
1041
1042         if (!qp1 || !qp2)
1043                 goto error;
1044
1045         if (qp1->div->n_row < qp2->div->n_row)
1046                 return isl_qpolynomial_add(qp2, qp1);
1047
1048         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1049         if (!compatible_divs(qp1->div, qp2->div))
1050                 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1051
1052         qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1053         if (!qp1->upoly)
1054                 goto error;
1055
1056         isl_qpolynomial_free(qp2);
1057
1058         return qp1;
1059 error:
1060         isl_qpolynomial_free(qp1);
1061         isl_qpolynomial_free(qp2);
1062         return NULL;
1063 }
1064
1065 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1066         __isl_take isl_qpolynomial *qp2)
1067 {
1068         return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1069 }
1070
1071 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1072 {
1073         qp = isl_qpolynomial_cow(qp);
1074
1075         if (!qp)
1076                 return NULL;
1077
1078         qp->upoly = isl_upoly_neg(qp->upoly);
1079         if (!qp->upoly)
1080                 goto error;
1081
1082         return qp;
1083 error:
1084         isl_qpolynomial_free(qp);
1085         return NULL;
1086 }
1087
1088 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1089         __isl_take isl_qpolynomial *qp2)
1090 {
1091         qp1 = isl_qpolynomial_cow(qp1);
1092
1093         if (!qp1 || !qp2)
1094                 goto error;
1095
1096         if (qp1->div->n_row < qp2->div->n_row)
1097                 return isl_qpolynomial_mul(qp2, qp1);
1098
1099         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1100         if (!compatible_divs(qp1->div, qp2->div))
1101                 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1102
1103         qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1104         if (!qp1->upoly)
1105                 goto error;
1106
1107         isl_qpolynomial_free(qp2);
1108
1109         return qp1;
1110 error:
1111         isl_qpolynomial_free(qp1);
1112         isl_qpolynomial_free(qp2);
1113         return NULL;
1114 }
1115
1116 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1117 {
1118         struct isl_qpolynomial *qp;
1119         struct isl_upoly_cst *cst;
1120
1121         qp = isl_qpolynomial_alloc(dim, 0);
1122         if (!qp)
1123                 return NULL;
1124
1125         qp->upoly = isl_upoly_zero(dim->ctx);
1126         if (!qp->upoly)
1127                 goto error;
1128
1129         return qp;
1130 error:
1131         isl_qpolynomial_free(qp);
1132         return NULL;
1133 }
1134
1135 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1136 {
1137         struct isl_qpolynomial *qp;
1138         struct isl_upoly_cst *cst;
1139
1140         qp = isl_qpolynomial_alloc(dim, 0);
1141         if (!qp)
1142                 return NULL;
1143
1144         qp->upoly = isl_upoly_infty(dim->ctx);
1145         if (!qp->upoly)
1146                 goto error;
1147
1148         return qp;
1149 error:
1150         isl_qpolynomial_free(qp);
1151         return NULL;
1152 }
1153
1154 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1155 {
1156         struct isl_qpolynomial *qp;
1157         struct isl_upoly_cst *cst;
1158
1159         qp = isl_qpolynomial_alloc(dim, 0);
1160         if (!qp)
1161                 return NULL;
1162
1163         qp->upoly = isl_upoly_nan(dim->ctx);
1164         if (!qp->upoly)
1165                 goto error;
1166
1167         return qp;
1168 error:
1169         isl_qpolynomial_free(qp);
1170         return NULL;
1171 }
1172
1173 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1174         isl_int v)
1175 {
1176         struct isl_qpolynomial *qp;
1177         struct isl_upoly_cst *cst;
1178
1179         qp = isl_qpolynomial_alloc(dim, 0);
1180         if (!qp)
1181                 return NULL;
1182
1183         qp->upoly = isl_upoly_zero(dim->ctx);
1184         if (!qp->upoly)
1185                 goto error;
1186         cst = isl_upoly_as_cst(qp->upoly);
1187         isl_int_set(cst->n, v);
1188
1189         return qp;
1190 error:
1191         isl_qpolynomial_free(qp);
1192         return NULL;
1193 }
1194
1195 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1196         isl_int *n, isl_int *d)
1197 {
1198         struct isl_upoly_cst *cst;
1199
1200         if (!qp)
1201                 return -1;
1202
1203         if (!isl_upoly_is_cst(qp->upoly))
1204                 return 0;
1205
1206         cst = isl_upoly_as_cst(qp->upoly);
1207         if (!cst)
1208                 return -1;
1209
1210         if (n)
1211                 isl_int_set(*n, cst->n);
1212         if (d)
1213                 isl_int_set(*d, cst->d);
1214
1215         return 1;
1216 }
1217
1218 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
1219         __isl_keep isl_qpolynomial *qp2)
1220 {
1221         if (!qp1 || !qp2)
1222                 return -1;
1223
1224         return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1225 }
1226
1227 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1228 {
1229         int i;
1230         struct isl_upoly_rec *rec;
1231
1232         if (isl_upoly_is_cst(up)) {
1233                 struct isl_upoly_cst *cst;
1234                 cst = isl_upoly_as_cst(up);
1235                 if (!cst)
1236                         return;
1237                 isl_int_lcm(*d, *d, cst->d);
1238                 return;
1239         }
1240
1241         rec = isl_upoly_as_rec(up);
1242         if (!rec)
1243                 return;
1244
1245         for (i = 0; i < rec->n; ++i)
1246                 upoly_update_den(rec->p[i], d);
1247 }
1248
1249 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1250 {
1251         isl_int_set_si(*d, 1);
1252         if (!qp)
1253                 return;
1254         upoly_update_den(qp->upoly, d);
1255 }
1256
1257 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1258         int pos, int power)
1259 {
1260         int i;
1261         struct isl_qpolynomial *qp;
1262         struct isl_upoly_rec *rec;
1263         struct isl_upoly_cst *cst;
1264         struct isl_ctx *ctx;
1265
1266         if (!dim)
1267                 return NULL;
1268
1269         ctx = dim->ctx;
1270
1271         qp = isl_qpolynomial_alloc(dim, 0);
1272         if (!qp)
1273                 return NULL;
1274
1275         qp->upoly = &isl_upoly_alloc_rec(ctx, pos, 1 + power)->up;
1276         if (!qp->upoly)
1277                 goto error;
1278         rec = isl_upoly_as_rec(qp->upoly);
1279         for (i = 0; i < 1 + power; ++i) {
1280                 rec->p[i] = isl_upoly_zero(ctx);
1281                 if (!rec->p[i])
1282                         goto error;
1283                 rec->n++;
1284         }
1285         cst = isl_upoly_as_cst(rec->p[power]);
1286         isl_int_set_si(cst->n, 1);
1287
1288         return qp;
1289 error:
1290         isl_qpolynomial_free(qp);
1291         return NULL;
1292 }
1293
1294 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1295         enum isl_dim_type type, unsigned pos)
1296 {
1297         if (!dim)
1298                 return NULL;
1299
1300         isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1301         isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1302
1303         if (type == isl_dim_set)
1304                 pos += isl_dim_size(dim, isl_dim_param);
1305
1306         return isl_qpolynomial_pow(dim, pos, 1);
1307 error:
1308         isl_dim_free(dim);
1309         return NULL;
1310 }
1311
1312 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1313         int power)
1314 {
1315         struct isl_qpolynomial *qp = NULL;
1316         struct isl_upoly_rec *rec;
1317         struct isl_upoly_cst *cst;
1318         int i;
1319         int pos;
1320
1321         if (!div)
1322                 return NULL;
1323         isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1324
1325         qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1);
1326         if (!qp)
1327                 goto error;
1328
1329         isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1330         isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1331
1332         pos = isl_dim_total(qp->dim);
1333         qp->upoly = &isl_upoly_alloc_rec(div->ctx, pos, 1 + power)->up;
1334         if (!qp->upoly)
1335                 goto error;
1336         rec = isl_upoly_as_rec(qp->upoly);
1337         for (i = 0; i < 1 + power; ++i) {
1338                 rec->p[i] = isl_upoly_zero(div->ctx);
1339                 if (!rec->p[i])
1340                         goto error;
1341                 rec->n++;
1342         }
1343         cst = isl_upoly_as_cst(rec->p[power]);
1344         isl_int_set_si(cst->n, 1);
1345
1346         isl_div_free(div);
1347
1348         return qp;
1349 error:
1350         isl_qpolynomial_free(qp);
1351         isl_div_free(div);
1352         return NULL;
1353 }
1354
1355 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1356 {
1357         return isl_qpolynomial_div_pow(div, 1);
1358 }
1359
1360 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1361         const isl_int n, const isl_int d)
1362 {
1363         struct isl_qpolynomial *qp;
1364         struct isl_upoly_cst *cst;
1365
1366         qp = isl_qpolynomial_alloc(dim, 0);
1367         if (!qp)
1368                 return NULL;
1369
1370         qp->upoly = isl_upoly_zero(dim->ctx);
1371         if (!qp->upoly)
1372                 goto error;
1373         cst = isl_upoly_as_cst(qp->upoly);
1374         isl_int_set(cst->n, n);
1375         isl_int_set(cst->d, d);
1376
1377         return qp;
1378 error:
1379         isl_qpolynomial_free(qp);
1380         return NULL;
1381 }
1382
1383 #undef PW
1384 #define PW isl_pw_qpolynomial
1385 #undef EL
1386 #define EL isl_qpolynomial
1387 #undef IS_ZERO
1388 #define IS_ZERO is_zero
1389 #undef FIELD
1390 #define FIELD qp
1391 #undef ADD
1392 #define ADD add
1393
1394 #include <isl_pw_templ.c>
1395
1396 #undef PW
1397 #define PW isl_pw_qpolynomial_fold
1398 #undef EL
1399 #define EL isl_qpolynomial_fold
1400 #undef IS_ZERO
1401 #define IS_ZERO is_empty
1402 #undef FIELD
1403 #define FIELD fold
1404 #undef ADD
1405 #define ADD fold
1406
1407 #include <isl_pw_templ.c>
1408
1409 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1410 {
1411         if (!pwqp)
1412                 return -1;
1413
1414         if (pwqp->n != -1)
1415                 return 0;
1416
1417         if (!isl_set_fast_is_universe(pwqp->p[0].set))
1418                 return 0;
1419
1420         return isl_qpolynomial_is_one(pwqp->p[0].qp);
1421 }
1422
1423 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1424         __isl_take isl_pw_qpolynomial *pwqp1,
1425         __isl_take isl_pw_qpolynomial *pwqp2)
1426 {
1427         int i, j, n;
1428         struct isl_pw_qpolynomial *res;
1429         isl_set *set;
1430
1431         if (!pwqp1 || !pwqp2)
1432                 goto error;
1433
1434         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1435                         goto error);
1436
1437         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1438                 isl_pw_qpolynomial_free(pwqp2);
1439                 return pwqp1;
1440         }
1441
1442         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1443                 isl_pw_qpolynomial_free(pwqp1);
1444                 return pwqp2;
1445         }
1446
1447         if (isl_pw_qpolynomial_is_one(pwqp1)) {
1448                 isl_pw_qpolynomial_free(pwqp1);
1449                 return pwqp2;
1450         }
1451
1452         if (isl_pw_qpolynomial_is_one(pwqp2)) {
1453                 isl_pw_qpolynomial_free(pwqp2);
1454                 return pwqp1;
1455         }
1456
1457         n = pwqp1->n * pwqp2->n;
1458         res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
1459
1460         for (i = 0; i < pwqp1->n; ++i) {
1461                 for (j = 0; j < pwqp2->n; ++j) {
1462                         struct isl_set *common;
1463                         struct isl_qpolynomial *prod;
1464                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1465                                                 isl_set_copy(pwqp2->p[j].set));
1466                         if (isl_set_fast_is_empty(common)) {
1467                                 isl_set_free(common);
1468                                 continue;
1469                         }
1470
1471                         prod = isl_qpolynomial_mul(
1472                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1473                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1474
1475                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
1476                 }
1477         }
1478
1479         isl_pw_qpolynomial_free(pwqp1);
1480         isl_pw_qpolynomial_free(pwqp2);
1481
1482         return res;
1483 error:
1484         isl_pw_qpolynomial_free(pwqp1);
1485         isl_pw_qpolynomial_free(pwqp2);
1486         return NULL;
1487 }
1488
1489 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1490         __isl_take isl_pw_qpolynomial *pwqp)
1491 {
1492         int i, j, n;
1493         struct isl_pw_qpolynomial *res;
1494         isl_set *set;
1495
1496         if (!pwqp)
1497                 return NULL;
1498
1499         if (isl_pw_qpolynomial_is_zero(pwqp))
1500                 return pwqp;
1501
1502         pwqp = isl_pw_qpolynomial_cow(pwqp);
1503
1504         for (i = 0; i < pwqp->n; ++i) {
1505                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1506                 if (!pwqp->p[i].qp)
1507                         goto error;
1508         }
1509
1510         return pwqp;
1511 error:
1512         isl_pw_qpolynomial_free(pwqp);
1513         return NULL;
1514 }
1515
1516 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1517         __isl_take isl_pw_qpolynomial *pwqp1,
1518         __isl_take isl_pw_qpolynomial *pwqp2)
1519 {
1520         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1521 }
1522
1523 __isl_give struct isl_upoly *isl_upoly_eval(
1524         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1525 {
1526         int i;
1527         struct isl_upoly_rec *rec;
1528         struct isl_upoly *res;
1529         struct isl_upoly *base;
1530
1531         if (isl_upoly_is_cst(up)) {
1532                 isl_vec_free(vec);
1533                 return up;
1534         }
1535
1536         rec = isl_upoly_as_rec(up);
1537         if (!rec)
1538                 goto error;
1539
1540         isl_assert(up->ctx, rec->n >= 1, goto error);
1541
1542         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1543
1544         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1545                                 isl_vec_copy(vec));
1546
1547         for (i = rec->n - 2; i >= 0; --i) {
1548                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1549                 res = isl_upoly_sum(res,
1550                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1551                                                             isl_vec_copy(vec)));
1552         }
1553
1554         isl_upoly_free(base);
1555         isl_upoly_free(up);
1556         isl_vec_free(vec);
1557         return res;
1558 error:
1559         isl_upoly_free(up);
1560         isl_vec_free(vec);
1561         return NULL;
1562 }
1563
1564 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1565         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1566 {
1567         isl_vec *ext;
1568         struct isl_upoly *up;
1569         isl_dim *dim;
1570
1571         if (!qp || !pnt)
1572                 goto error;
1573         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1574
1575         if (qp->div->n_row == 0)
1576                 ext = isl_vec_copy(pnt->vec);
1577         else {
1578                 int i;
1579                 unsigned dim = isl_dim_total(qp->dim);
1580                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1581                 if (!ext)
1582                         goto error;
1583
1584                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1585                 for (i = 0; i < qp->div->n_row; ++i) {
1586                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1587                                                 1 + dim + i, &ext->el[1+dim+i]);
1588                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1589                                         qp->div->row[i][0]);
1590                 }
1591         }
1592
1593         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1594         if (!up)
1595                 goto error;
1596
1597         dim = isl_dim_copy(qp->dim);
1598         isl_qpolynomial_free(qp);
1599         isl_point_free(pnt);
1600
1601         qp = isl_qpolynomial_alloc(dim, 0);
1602         if (!qp)
1603                 isl_upoly_free(up);
1604         else
1605                 qp->upoly = up;
1606
1607         return qp;
1608 error:
1609         isl_qpolynomial_free(qp);
1610         isl_point_free(pnt);
1611         return NULL;
1612 }
1613
1614 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
1615         __isl_keep struct isl_upoly_cst *cst2)
1616 {
1617         int cmp;
1618         isl_int t;
1619         isl_int_init(t);
1620         isl_int_mul(t, cst1->n, cst2->d);
1621         isl_int_submul(t, cst2->n, cst1->d);
1622         cmp = isl_int_sgn(t);
1623         isl_int_clear(t);
1624         return cmp;
1625 }
1626
1627 static __isl_give isl_qpolynomial *qpolynomial_min(
1628         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1629 {
1630         struct isl_upoly_cst *cst1, *cst2;
1631         int cmp;
1632
1633         if (!qp1 || !qp2)
1634                 goto error;
1635         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1636         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1637         cst1 = isl_upoly_as_cst(qp1->upoly);
1638         cst2 = isl_upoly_as_cst(qp2->upoly);
1639         cmp = isl_upoly_cmp(cst1, cst2);
1640
1641         if (cmp <= 0) {
1642                 isl_qpolynomial_free(qp2);
1643         } else {
1644                 isl_qpolynomial_free(qp1);
1645                 qp1 = qp2;
1646         }
1647         return qp1;
1648 error:
1649         isl_qpolynomial_free(qp1);
1650         isl_qpolynomial_free(qp2);
1651         return NULL;
1652 }
1653
1654 static __isl_give isl_qpolynomial *qpolynomial_max(
1655         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1656 {
1657         struct isl_upoly_cst *cst1, *cst2;
1658         int cmp;
1659
1660         if (!qp1 || !qp2)
1661                 goto error;
1662         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1663         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1664         cst1 = isl_upoly_as_cst(qp1->upoly);
1665         cst2 = isl_upoly_as_cst(qp2->upoly);
1666         cmp = isl_upoly_cmp(cst1, cst2);
1667
1668         if (cmp >= 0) {
1669                 isl_qpolynomial_free(qp2);
1670         } else {
1671                 isl_qpolynomial_free(qp1);
1672                 qp1 = qp2;
1673         }
1674         return qp1;
1675 error:
1676         isl_qpolynomial_free(qp1);
1677         isl_qpolynomial_free(qp2);
1678         return NULL;
1679 }
1680
1681 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
1682         __isl_keep isl_qpolynomial_fold *fold2)
1683 {
1684         int i;
1685
1686         if (!fold1 || !fold2)
1687                 return -1;
1688
1689         if (fold1->n != fold2->n)
1690                 return 0;
1691
1692         /* We probably want to sort the qps first... */
1693         for (i = 0; i < fold1->n; ++i) {
1694                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
1695                 if (eq < 0 || !eq)
1696                         return eq;
1697         }
1698
1699         return 1;
1700 }
1701
1702 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
1703         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
1704 {
1705         isl_qpolynomial *qp;
1706
1707         if (!fold || !pnt)
1708                 goto error;
1709         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
1710         isl_assert(pnt->dim->ctx,
1711                 fold->type == isl_fold_max || fold->type == isl_fold_min,
1712                 goto error);
1713
1714         if (fold->n == 0)
1715                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
1716         else {
1717                 int i;
1718                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
1719                                                 isl_point_copy(pnt));
1720                 for (i = 1; i < fold->n; ++i) {
1721                         isl_qpolynomial *qp_i;
1722                         qp_i = isl_qpolynomial_eval(
1723                                             isl_qpolynomial_copy(fold->qp[i]),
1724                                             isl_point_copy(pnt));
1725                         if (fold->type == isl_fold_max)
1726                                 qp = qpolynomial_max(qp, qp_i);
1727                         else
1728                                 qp = qpolynomial_min(qp, qp_i);
1729                 }
1730         }
1731         isl_qpolynomial_fold_free(fold);
1732         isl_point_free(pnt);
1733
1734         return qp;
1735 error:
1736         isl_qpolynomial_fold_free(fold);
1737         isl_point_free(pnt);
1738         return NULL;
1739 }
1740
1741 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1742         enum isl_dim_type dst_type, unsigned dst_pos,
1743         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1744 {
1745         if (!qp)
1746                 return NULL;
1747
1748         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1749         if (!qp->dim)
1750                 goto error;
1751
1752         /* Do something to polynomials when needed; later */
1753
1754         return qp;
1755 error:
1756         isl_qpolynomial_free(qp);
1757         return NULL;
1758 }
1759
1760 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1761         __isl_take isl_pw_qpolynomial *pwqp,
1762         enum isl_dim_type dst_type, unsigned dst_pos,
1763         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1764 {
1765         int i;
1766
1767         if (!pwqp)
1768                 return NULL;
1769
1770         pwqp->dim = isl_dim_move(pwqp->dim,
1771                                     dst_type, dst_pos, src_type, src_pos, n);
1772         if (!pwqp->dim)
1773                 goto error;
1774
1775         for (i = 0; i < pwqp->n; ++i) {
1776                 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1777                                                 src_type, src_pos, n);
1778                 if (!pwqp->p[i].set)
1779                         goto error;
1780                 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1781                                         dst_type, dst_pos, src_type, src_pos, n);
1782                 if (!pwqp->p[i].qp)
1783                         goto error;
1784         }
1785
1786         return pwqp;
1787 error:
1788         isl_pw_qpolynomial_free(pwqp);
1789         return NULL;
1790 }
1791
1792 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
1793         __isl_keep isl_pw_qpolynomial *pwqp)
1794 {
1795         if (!pwqp)
1796                 return NULL;
1797
1798         return isl_dim_copy(pwqp->dim);
1799 }
1800
1801 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
1802         enum isl_dim_type type)
1803 {
1804         return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
1805 }
1806
1807 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
1808         __isl_take isl_mat *div)
1809 {
1810         isl_term *term;
1811         int n;
1812
1813         if (!dim || !div)
1814                 goto error;
1815
1816         n = isl_dim_total(dim) + div->n_row;
1817
1818         term = isl_calloc(dim->ctx, struct isl_term,
1819                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
1820         if (!term)
1821                 goto error;
1822
1823         term->ref = 1;
1824         term->dim = dim;
1825         term->div = div;
1826         isl_int_init(term->n);
1827         isl_int_init(term->d);
1828
1829         return term;
1830 error:
1831         isl_dim_free(dim);
1832         isl_mat_free(div);
1833         return NULL;
1834 }
1835
1836 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
1837 {
1838         if (!term)
1839                 return NULL;
1840
1841         term->ref++;
1842         return term;
1843 }
1844
1845 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
1846 {
1847         int i;
1848         isl_term *dup;
1849         unsigned total;
1850
1851         if (term)
1852                 return NULL;
1853
1854         total = isl_dim_total(term->dim) + term->div->n_row;
1855
1856         dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
1857         if (!dup)
1858                 return NULL;
1859
1860         isl_int_set(dup->n, term->n);
1861         isl_int_set(dup->d, term->d);
1862
1863         for (i = 0; i < total; ++i)
1864                 dup->pow[i] = term->pow[i];
1865
1866         return dup;
1867 }
1868
1869 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
1870 {
1871         if (!term)
1872                 return NULL;
1873
1874         if (term->ref == 1)
1875                 return term;
1876         term->ref--;
1877         return isl_term_dup(term);
1878 }
1879
1880 void isl_term_free(__isl_take isl_term *term)
1881 {
1882         if (!term)
1883                 return;
1884
1885         if (--term->ref > 0)
1886                 return;
1887
1888         isl_dim_free(term->dim);
1889         isl_mat_free(term->div);
1890         isl_int_clear(term->n);
1891         isl_int_clear(term->d);
1892         free(term);
1893 }
1894
1895 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
1896 {
1897         if (!term)
1898                 return 0;
1899
1900         switch (type) {
1901         case isl_dim_param:
1902         case isl_dim_in:
1903         case isl_dim_out:       return isl_dim_size(term->dim, type);
1904         case isl_dim_div:       return term->div->n_row;
1905         case isl_dim_all:       return isl_dim_total(term->dim) + term->div->n_row;
1906         }
1907 }
1908
1909 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
1910 {
1911         return term ? term->dim->ctx : NULL;
1912 }
1913
1914 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
1915 {
1916         if (!term)
1917                 return;
1918         isl_int_set(*n, term->n);
1919 }
1920
1921 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
1922 {
1923         if (!term)
1924                 return;
1925         isl_int_set(*d, term->d);
1926 }
1927
1928 int isl_term_get_exp(__isl_keep isl_term *term,
1929         enum isl_dim_type type, unsigned pos)
1930 {
1931         if (!term)
1932                 return -1;
1933
1934         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
1935
1936         if (type >= isl_dim_set)
1937                 pos += isl_dim_size(term->dim, isl_dim_param);
1938         if (type >= isl_dim_div)
1939                 pos += isl_dim_size(term->dim, isl_dim_set);
1940
1941         return term->pow[pos];
1942 }
1943
1944 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
1945 {
1946         isl_basic_map *bmap;
1947         unsigned total;
1948         int k;
1949
1950         if (!term)
1951                 return NULL;
1952
1953         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
1954                         return NULL);
1955
1956         total = term->div->n_col - term->div->n_row - 2;
1957         /* No nested divs for now */
1958         isl_assert(term->dim->ctx,
1959                 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
1960                                         term->div->n_row) == -1,
1961                 return NULL);
1962
1963         bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
1964         if ((k = isl_basic_map_alloc_div(bmap)) < 0)
1965                 goto error;
1966
1967         isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
1968
1969         return isl_basic_map_div(bmap, k);
1970 error:
1971         isl_basic_map_free(bmap);
1972         return NULL;
1973 }
1974
1975 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
1976         int (*fn)(__isl_take isl_term *term, void *user),
1977         __isl_take isl_term *term, void *user)
1978 {
1979         int i;
1980         struct isl_upoly_rec *rec;
1981
1982         if (!up || !term)
1983                 goto error;
1984
1985         if (isl_upoly_is_zero(up))
1986                 return term;
1987
1988         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
1989         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
1990         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
1991
1992         if (isl_upoly_is_cst(up)) {
1993                 struct isl_upoly_cst *cst;
1994                 cst = isl_upoly_as_cst(up);
1995                 if (!cst)
1996                         goto error;
1997                 term = isl_term_cow(term);
1998                 if (!term)
1999                         goto error;
2000                 isl_int_set(term->n, cst->n);
2001                 isl_int_set(term->d, cst->d);
2002                 if (fn(isl_term_copy(term), user) < 0)
2003                         goto error;
2004                 return term;
2005         }
2006
2007         rec = isl_upoly_as_rec(up);
2008         if (!rec)
2009                 goto error;
2010
2011         for (i = 0; i < rec->n; ++i) {
2012                 term = isl_term_cow(term);
2013                 if (!term)
2014                         goto error;
2015                 term->pow[up->var] = i;
2016                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2017                 if (!term)
2018                         goto error;
2019         }
2020         term->pow[up->var] = 0;
2021
2022         return term;
2023 error:
2024         isl_term_free(term);
2025         return NULL;
2026 }
2027
2028 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
2029         int (*fn)(__isl_take isl_term *term, void *user), void *user)
2030 {
2031         isl_term *term;
2032
2033         if (!qp)
2034                 return -1;
2035
2036         term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
2037         if (!term)
2038                 return -1;
2039
2040         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
2041
2042         isl_term_free(term);
2043
2044         return term ? 0 : -1;
2045 }
2046
2047 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2048         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2049                     void *user), void *user)
2050 {
2051         int i;
2052
2053         if (!pwqp)
2054                 return -1;
2055
2056         for (i = 0; i < pwqp->n; ++i)
2057                 if (fn(isl_set_copy(pwqp->p[i].set),
2058                                 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
2059                         return -1;
2060
2061         return 0;
2062 }
2063
2064 static int any_divs(__isl_keep isl_set *set)
2065 {
2066         int i;
2067
2068         if (!set)
2069                 return -1;
2070
2071         for (i = 0; i < set->n; ++i)
2072                 if (set->p[i]->n_div > 0)
2073                         return 1;
2074
2075         return 0;
2076 }
2077
2078 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
2079         __isl_take isl_dim *dim)
2080 {
2081         if (!qp || !dim)
2082                 goto error;
2083
2084         if (isl_dim_equal(qp->dim, dim)) {
2085                 isl_dim_free(dim);
2086                 return qp;
2087         }
2088
2089         qp = isl_qpolynomial_cow(qp);
2090         if (!qp)
2091                 goto error;
2092
2093         if (qp->div->n_row) {
2094                 int i;
2095                 int extra;
2096                 unsigned total;
2097                 int *exp;
2098
2099                 extra = isl_dim_size(dim, isl_dim_set) -
2100                                 isl_dim_size(qp->dim, isl_dim_set);
2101                 total = isl_dim_total(qp->dim);
2102                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
2103                 if (!exp)
2104                         goto error;
2105                 for (i = 0; i < qp->div->n_row; ++i)
2106                         exp[i] = extra + i;
2107                 qp->upoly = expand(qp->upoly, exp, total);
2108                 free(exp);
2109                 if (!qp->upoly)
2110                         goto error;
2111                 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
2112                 if (!qp->div)
2113                         goto error;
2114                 for (i = 0; i < qp->div->n_row; ++i)
2115                         isl_seq_clr(qp->div->row[i] + 2 + total, extra);
2116         }
2117
2118         isl_dim_free(qp->dim);
2119         qp->dim = dim;
2120
2121         return qp;
2122 error:
2123         isl_dim_free(dim);
2124         isl_qpolynomial_free(qp);
2125         return NULL;
2126 }
2127
2128 static int foreach_lifted_subset(__isl_take isl_set *set,
2129         __isl_take isl_qpolynomial *qp,
2130         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2131                     void *user), void *user)
2132 {
2133         int i;
2134
2135         if (!set || !qp)
2136                 goto error;
2137
2138         for (i = 0; i < set->n; ++i) {
2139                 isl_set *lift;
2140                 isl_qpolynomial *copy;
2141
2142                 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2143                 lift = isl_set_lift(lift);
2144
2145                 copy = isl_qpolynomial_copy(qp);
2146                 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2147
2148                 if (fn(lift, copy, user) < 0)
2149                         goto error;
2150         }
2151
2152         isl_set_free(set);
2153         isl_qpolynomial_free(qp);
2154
2155         return 0;
2156 error:
2157         isl_set_free(set);
2158         isl_qpolynomial_free(qp);
2159         return -1;
2160 }
2161
2162 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2163         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2164                     void *user), void *user)
2165 {
2166         int i;
2167
2168         if (!pwqp)
2169                 return -1;
2170
2171         for (i = 0; i < pwqp->n; ++i) {
2172                 isl_set *set;
2173                 isl_qpolynomial *qp;
2174
2175                 set = isl_set_copy(pwqp->p[i].set);
2176                 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2177                 if (!any_divs(set)) {
2178                         if (fn(set, qp, user) < 0)
2179                                 return -1;
2180                         continue;
2181                 }
2182                 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2183                         return -1;
2184         }
2185
2186         return 0;
2187 }
2188
2189 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
2190         enum isl_fold type, __isl_take isl_dim *dim, int n)
2191 {
2192         isl_qpolynomial_fold *fold;
2193
2194         if (!dim)
2195                 goto error;
2196
2197         isl_assert(dim->ctx, n >= 0, goto error);
2198         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
2199                         sizeof(struct isl_qpolynomial_fold) +
2200                         (n - 1) * sizeof(struct isl_qpolynomial *));
2201         if (!fold)
2202                 goto error;
2203
2204         fold->ref = 1;
2205         fold->size = n;
2206         fold->n = 0;
2207         fold->type = type;
2208         fold->dim = dim;
2209
2210         return fold;
2211 error:
2212         isl_dim_free(dim);
2213         return NULL;
2214 }
2215
2216 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
2217         __isl_take isl_dim *dim)
2218 {
2219         return qpolynomial_fold_alloc(type, dim, 0);
2220 }
2221
2222 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
2223         enum isl_fold type, __isl_take isl_qpolynomial *qp)
2224 {
2225         isl_qpolynomial_fold *fold;
2226
2227         if (!qp)
2228                 return NULL;
2229
2230         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
2231         if (!fold)
2232                 goto error;
2233
2234         fold->qp[0] = qp;
2235         fold->n++;
2236
2237         return fold;
2238 error:
2239         isl_qpolynomial_fold_free(fold);
2240         isl_qpolynomial_free(qp);
2241         return NULL;
2242 }
2243
2244 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
2245         __isl_keep isl_qpolynomial_fold *fold)
2246 {
2247         if (!fold)
2248                 return NULL;
2249
2250         fold->ref++;
2251         return fold;
2252 }
2253
2254 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
2255 {
2256         int i;
2257
2258         if (!fold)
2259                 return;
2260         if (--fold->ref > 0)
2261                 return;
2262
2263         for (i = 0; i < fold->n; ++i)
2264                 isl_qpolynomial_free(fold->qp[i]);
2265         isl_dim_free(fold->dim);
2266         free(fold);
2267 }
2268
2269 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
2270 {
2271         if (!fold)
2272                 return -1;
2273
2274         return fold->n == 0;
2275 }
2276
2277 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
2278         __isl_take isl_qpolynomial_fold *fold1,
2279         __isl_take isl_qpolynomial_fold *fold2)
2280 {
2281         int i;
2282         struct isl_qpolynomial_fold *res = NULL;
2283
2284         if (!fold1 || !fold2)
2285                 goto error;
2286
2287         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
2288         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
2289                         goto error);
2290
2291         if (isl_qpolynomial_fold_is_empty(fold1)) {
2292                 isl_qpolynomial_fold_free(fold1);
2293                 return fold2;
2294         }
2295
2296         if (isl_qpolynomial_fold_is_empty(fold2)) {
2297                 isl_qpolynomial_fold_free(fold2);
2298                 return fold1;
2299         }
2300
2301         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
2302                                         fold1->n + fold2->n);
2303         if (!res)
2304                 goto error;
2305
2306         for (i = 0; i < fold1->n; ++i) {
2307                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
2308                 if (!res->qp[res->n])
2309                         goto error;
2310                 res->n++;
2311         }
2312
2313         for (i = 0; i < fold2->n; ++i) {
2314                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
2315                 if (!res->qp[res->n])
2316                         goto error;
2317                 res->n++;
2318         }
2319
2320         isl_qpolynomial_fold_free(fold1);
2321         isl_qpolynomial_fold_free(fold2);
2322
2323         return res;
2324 error:
2325         isl_qpolynomial_fold_free(res);
2326         isl_qpolynomial_fold_free(fold1);
2327         isl_qpolynomial_fold_free(fold2);
2328         return NULL;
2329 }
2330
2331 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2332         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
2333 {
2334         int i;
2335         isl_pw_qpolynomial_fold *pwf;
2336
2337         if (!pwqp)
2338                 return NULL;
2339
2340         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
2341
2342         for (i = 0; i < pwqp->n; ++i)
2343                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
2344                         isl_set_copy(pwqp->p[i].set),
2345                         isl_qpolynomial_fold_alloc(type,
2346                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
2347
2348         isl_pw_qpolynomial_free(pwqp);
2349
2350         return pwf;
2351 }
2352
2353 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2354 {
2355         struct isl_upoly_rec *rec;
2356         int i;
2357
2358         if (!up)
2359                 return -1;
2360
2361         if (isl_upoly_is_cst(up))
2362                 return 0;
2363
2364         if (up->var < d)
2365                 active[up->var] = 1;
2366
2367         rec = isl_upoly_as_rec(up);
2368         for (i = 0; i < rec->n; ++i)
2369                 if (up_set_active(rec->p[i], active, d) < 0)
2370                         return -1;
2371
2372         return 0;
2373 }
2374
2375 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2376 {
2377         int i, j;
2378         int d = isl_dim_total(qp->dim);
2379
2380         for (i = 0; i < d; ++i)
2381                 for (j = 0; j < qp->div->n_row; ++j) {
2382                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
2383                                 continue;
2384                         active[i] = 1;
2385                         break;
2386                 }
2387
2388         return up_set_active(qp->upoly, active, d);
2389 }
2390
2391 /* For each parameter or variable that does not appear in qp,
2392  * first eliminate the variable from all constraints and then set it to zero.
2393  */
2394 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
2395         __isl_keep isl_qpolynomial *qp)
2396 {
2397         int *active = NULL;
2398         int i;
2399         int d;
2400         unsigned nparam;
2401         unsigned nvar;
2402
2403         if (!set || !qp)
2404                 goto error;
2405
2406         d = isl_dim_total(set->dim);
2407         active = isl_calloc_array(set->ctx, int, d);
2408         if (set_active(qp, active) < 0)
2409                 goto error;
2410
2411         for (i = 0; i < d; ++i)
2412                 if (!active[i])
2413                         break;
2414
2415         if (i == d) {
2416                 free(active);
2417                 return set;
2418         }
2419
2420         nparam = isl_dim_size(set->dim, isl_dim_param);
2421         nvar = isl_dim_size(set->dim, isl_dim_set);
2422         for (i = 0; i < nparam; ++i) {
2423                 if (active[i])
2424                         continue;
2425                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
2426                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
2427         }
2428         for (i = 0; i < nvar; ++i) {
2429                 if (active[i])
2430                         continue;
2431                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
2432                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
2433         }
2434
2435         free(active);
2436
2437         return set;
2438 error:
2439         free(active);
2440         isl_set_free(set);
2441         return NULL;
2442 }
2443
2444 struct isl_max_data {
2445         isl_qpolynomial *qp;
2446         int first;
2447         isl_qpolynomial *max;
2448 };
2449
2450 static int max_fn(__isl_take isl_point *pnt, void *user)
2451 {
2452         struct isl_max_data *data = (struct isl_max_data *)user;
2453         isl_qpolynomial *val;
2454
2455         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
2456         if (data->first) {
2457                 data->first = 0;
2458                 data->max = val;
2459         } else {
2460                 data->max = qpolynomial_max(data->max, val);
2461         }
2462
2463         return 0;
2464 }
2465
2466 static __isl_give isl_qpolynomial *guarded_qpolynomial_max(
2467         __isl_take isl_set *set, __isl_take isl_qpolynomial *qp)
2468 {
2469         struct isl_max_data data = { NULL, 1, NULL };
2470
2471         if (!set || !qp)
2472                 goto error;
2473
2474         if (isl_upoly_is_cst(qp->upoly)) {
2475                 isl_set_free(set);
2476                 return qp;
2477         }
2478
2479         set = fix_inactive(set, qp);
2480
2481         data.qp = qp;
2482         if (isl_set_foreach_point(set, max_fn, &data) < 0)
2483                 goto error;
2484
2485         isl_set_free(set);
2486         isl_qpolynomial_free(qp);
2487         return data.max;
2488 error:
2489         isl_set_free(set);
2490         isl_qpolynomial_free(qp);
2491         isl_qpolynomial_free(data.max);
2492         return NULL;
2493 }
2494
2495 /* Compute the maximal value attained by the piecewise quasipolynomial
2496  * on its domain or zero if the domain is empty.
2497  * In the worst case, the domain is scanned completely,
2498  * so the domain is assumed to be bounded.
2499  */
2500 __isl_give isl_qpolynomial *isl_pw_qpolynomial_max(
2501         __isl_take isl_pw_qpolynomial *pwqp)
2502 {
2503         int i;
2504         isl_qpolynomial *max;
2505
2506         if (!pwqp)
2507                 return NULL;
2508
2509         if (pwqp->n == 0) {
2510                 isl_dim *dim = isl_dim_copy(pwqp->dim);
2511                 isl_pw_qpolynomial_free(pwqp);
2512                 return isl_qpolynomial_zero(dim);
2513         }
2514
2515         max = guarded_qpolynomial_max(isl_set_copy(pwqp->p[0].set),
2516                                         isl_qpolynomial_copy(pwqp->p[0].qp));
2517         for (i = 1; i < pwqp->n; ++i) {
2518                 isl_qpolynomial *max_i;
2519                 max_i = guarded_qpolynomial_max(isl_set_copy(pwqp->p[i].set),
2520                                             isl_qpolynomial_copy(pwqp->p[i].qp));
2521                 max = qpolynomial_max(max, max_i);
2522         }
2523
2524         isl_pw_qpolynomial_free(pwqp);
2525         return max;
2526 }