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         cst = isl_upoly_copy(rec->p[0]);
 471         isl_upoly_free(up);
 472         return cst;
 473 error:
 474         isl_upoly_free(up);
 475         return NULL;
 476 }
 477
 478 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
 479         __isl_take struct isl_upoly *up2)
 480 {
 481         int i;
 482         struct isl_upoly_rec *rec1, *rec2;
 483
 484         if (!up1 || !up2)
 485                 goto error;
 486
 487         if (isl_upoly_is_nan(up1)) {
 488                 isl_upoly_free(up2);
 489                 return up1;
 490         }
 491
 492         if (isl_upoly_is_nan(up2)) {
 493                 isl_upoly_free(up1);
 494                 return up2;
 495         }
 496
 497         if (isl_upoly_is_zero(up1)) {
 498                 isl_upoly_free(up1);
 499                 return up2;
 500         }
 501
 502         if (isl_upoly_is_zero(up2)) {
 503                 isl_upoly_free(up2);
 504                 return up1;
 505         }
 506
 507         if (up1->var < up2->var)
 508                 return isl_upoly_sum(up2, up1);
 509
 510         if (up2->var < up1->var) {
 511                 struct isl_upoly_rec *rec;
 512                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 513                         isl_upoly_free(up1);
 514                         return up2;
 515                 }
 516                 up1 = isl_upoly_cow(up1);
 517                 rec = isl_upoly_as_rec(up1);
 518                 if (!rec)
 519                         goto error;
 520                 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
 521                 if (rec->n == 1)
 522                         up1 = replace_by_constant_term(up1);
 523                 return up1;
 524         }
 525
 526         if (isl_upoly_is_cst(up1))
 527                 return isl_upoly_sum_cst(up1, up2);
 528
 529         rec1 = isl_upoly_as_rec(up1);
 530         rec2 = isl_upoly_as_rec(up2);
 531         if (!rec1 || !rec2)
 532                 goto error;
 533
 534         if (rec1->n < rec2->n)
 535                 return isl_upoly_sum(up2, up1);
 536
 537         up1 = isl_upoly_cow(up1);
 538         rec1 = isl_upoly_as_rec(up1);
 539         if (!rec1)
 540                 goto error;
 541
 542         for (i = rec2->n - 1; i >= 0; --i) {
 543                 rec1->p[i] = isl_upoly_sum(rec1->p[i],
 544                                             isl_upoly_copy(rec2->p[i]));
 545                 if (!rec1->p[i])
 546                         goto error;
 547                 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
 548                         isl_upoly_free(rec1->p[i]);
 549                         rec1->n--;
 550                 }
 551         }
 552
 553         if (rec1->n == 0)
 554                 up1 = replace_by_zero(up1);
 555         else if (rec1->n == 1)
 556                 up1 = replace_by_constant_term(up1);
 557
 558         isl_upoly_free(up2);
 559
 560         return up1;
 561 error:
 562         isl_upoly_free(up1);
 563         isl_upoly_free(up2);
 564         return NULL;
 565 }
 566
 567 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
 568 {
 569         struct isl_upoly_cst *cst;
 570
 571         if (isl_upoly_is_zero(up))
 572                 return up;
 573
 574         up = isl_upoly_cow(up);
 575         if (!up)
 576                 return NULL;
 577
 578         cst = isl_upoly_as_cst(up);
 579
 580         isl_int_neg(cst->n, cst->n);
 581
 582         return up;
 583 }
 584
 585 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
 586 {
 587         int i;
 588         struct isl_upoly_rec *rec;
 589
 590         if (!up)
 591                 return NULL;
 592
 593         if (isl_upoly_is_cst(up))
 594                 return isl_upoly_neg_cst(up);
 595
 596         up = isl_upoly_cow(up);
 597         rec = isl_upoly_as_rec(up);
 598         if (!rec)
 599                 goto error;
 600
 601         for (i = 0; i < rec->n; ++i) {
 602                 rec->p[i] = isl_upoly_neg(rec->p[i]);
 603                 if (!rec->p[i])
 604                         goto error;
 605         }
 606
 607         return up;
 608 error:
 609         isl_upoly_free(up);
 610         return NULL;
 611 }
 612
 613 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
 614         __isl_take struct isl_upoly *up2)
 615 {
 616         struct isl_upoly_cst *cst1;
 617         struct isl_upoly_cst *cst2;
 618
 619         up1 = isl_upoly_cow(up1);
 620         if (!up1 || !up2)
 621                 goto error;
 622
 623         cst1 = isl_upoly_as_cst(up1);
 624         cst2 = isl_upoly_as_cst(up2);
 625
 626         isl_int_mul(cst1->n, cst1->n, cst2->n);
 627         isl_int_mul(cst1->d, cst1->d, cst2->d);
 628
 629         isl_upoly_cst_reduce(cst1);
 630
 631         isl_upoly_free(up2);
 632         return up1;
 633 error:
 634         isl_upoly_free(up1);
 635         isl_upoly_free(up2);
 636         return NULL;
 637 }
 638
 639 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
 640         __isl_take struct isl_upoly *up2)
 641 {
 642         struct isl_upoly_rec *rec1;
 643         struct isl_upoly_rec *rec2;
 644         struct isl_upoly_rec *res;
 645         int i, j;
 646         int size;
 647
 648         rec1 = isl_upoly_as_rec(up1);
 649         rec2 = isl_upoly_as_rec(up2);
 650         if (!rec1 || !rec2)
 651                 goto error;
 652         size = rec1->n + rec2->n - 1;
 653         res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
 654         if (!res)
 655                 goto error;
 656
 657         for (i = 0; i < rec1->n; ++i) {
 658                 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
 659                                             isl_upoly_copy(rec1->p[i]));
 660                 if (!res->p[i])
 661                         goto error;
 662                 res->n++;
 663         }
 664         for (; i < size; ++i) {
 665                 res->p[i] = isl_upoly_zero(up1->ctx);
 666                 if (!res->p[i])
 667                         goto error;
 668                 res->n++;
 669         }
 670         for (i = 0; i < rec1->n; ++i) {
 671                 for (j = 1; j < rec2->n; ++j) {
 672                         struct isl_upoly *up;
 673                         up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
 674                                             isl_upoly_copy(rec1->p[i]));
 675                         res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
 676                         if (!res->p[i + j])
 677                                 goto error;
 678                 }
 679         }
 680
 681         isl_upoly_free(up1);
 682         isl_upoly_free(up2);
 683
 684         return &res->up;
 685 error:
 686         isl_upoly_free(up1);
 687         isl_upoly_free(up2);
 688         isl_upoly_free(&res->up);
 689         return NULL;
 690 }
 691
 692 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
 693         __isl_take struct isl_upoly *up2)
 694 {
 695         if (!up1 || !up2)
 696                 goto error;
 697
 698         if (isl_upoly_is_nan(up1)) {
 699                 isl_upoly_free(up2);
 700                 return up1;
 701         }
 702
 703         if (isl_upoly_is_nan(up2)) {
 704                 isl_upoly_free(up1);
 705                 return up2;
 706         }
 707
 708         if (isl_upoly_is_zero(up1)) {
 709                 isl_upoly_free(up2);
 710                 return up1;
 711         }
 712
 713         if (isl_upoly_is_zero(up2)) {
 714                 isl_upoly_free(up1);
 715                 return up2;
 716         }
 717
 718         if (isl_upoly_is_one(up1)) {
 719                 isl_upoly_free(up1);
 720                 return up2;
 721         }
 722
 723         if (isl_upoly_is_one(up2)) {
 724                 isl_upoly_free(up2);
 725                 return up1;
 726         }
 727
 728         if (up1->var < up2->var)
 729                 return isl_upoly_mul(up2, up1);
 730
 731         if (up2->var < up1->var) {
 732                 int i;
 733                 struct isl_upoly_rec *rec;
 734                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 735                         isl_ctx *ctx = up1->ctx;
 736                         isl_upoly_free(up1);
 737                         isl_upoly_free(up2);
 738                         return isl_upoly_nan(ctx);
 739                 }
 740                 up1 = isl_upoly_cow(up1);
 741                 rec = isl_upoly_as_rec(up1);
 742                 if (!rec)
 743                         goto error;
 744
 745                 for (i = 0; i < rec->n; ++i) {
 746                         rec->p[i] = isl_upoly_mul(rec->p[i],
 747                                                     isl_upoly_copy(up2));
 748                         if (!rec->p[i])
 749                                 goto error;
 750                 }
 751                 isl_upoly_free(up2);
 752                 return up1;
 753         }
 754
 755         if (isl_upoly_is_cst(up1))
 756                 return isl_upoly_mul_cst(up1, up2);
 757
 758         return isl_upoly_mul_rec(up1, up2);
 759 error:
 760         isl_upoly_free(up1);
 761         isl_upoly_free(up2);
 762         return NULL;
 763 }
 764
 765 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
 766         unsigned n_div)
 767 {
 768         struct isl_qpolynomial *qp;
 769         unsigned total;
 770
 771         if (!dim)
 772                 return NULL;
 773
 774         total = isl_dim_total(dim);
 775
 776         qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
 777         if (!qp)
 778                 return NULL;
 779
 780         qp->ref = 1;
 781         qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
 782         if (!qp->div)
 783                 goto error;
 784
 785         qp->dim = dim;
 786
 787         return qp;
 788 error:
 789         isl_dim_free(dim);
 790         isl_qpolynomial_free(qp);
 791         return NULL;
 792 }
 793
 794 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
 795 {
 796         if (!qp)
 797                 return NULL;
 798
 799         qp->ref++;
 800         return qp;
 801 }
 802
 803 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
 804 {
 805         struct isl_qpolynomial *dup;
 806
 807         if (!qp)
 808                 return NULL;
 809
 810         dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row);
 811         if (!dup)
 812                 return NULL;
 813         isl_mat_free(dup->div);
 814         dup->div = isl_mat_copy(qp->div);
 815         if (!dup->div)
 816                 goto error;
 817         dup->upoly = isl_upoly_copy(qp->upoly);
 818         if (!dup->upoly)
 819                 goto error;
 820
 821         return dup;
 822 error:
 823         isl_qpolynomial_free(dup);
 824         return NULL;
 825 }
 826
 827 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
 828 {
 829         if (!qp)
 830                 return NULL;
 831
 832         if (qp->ref == 1)
 833                 return qp;
 834         qp->ref--;
 835         return isl_qpolynomial_dup(qp);
 836 }
 837
 838 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
 839 {
 840         if (!qp)
 841                 return;
 842
 843         if (--qp->ref > 0)
 844                 return;
 845
 846         isl_dim_free(qp->dim);
 847         isl_mat_free(qp->div);
 848         isl_upoly_free(qp->upoly);
 849
 850         free(qp);
 851 }
 852
 853 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
 854 {
 855         int n_row, n_col;
 856         int equal;
 857
 858         isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
 859                                 div1->n_col >= div2->n_col, return -1);
 860
 861         if (div1->n_row == div2->n_row)
 862                 return isl_mat_is_equal(div1, div2);
 863
 864         n_row = div1->n_row;
 865         n_col = div1->n_col;
 866         div1->n_row = div2->n_row;
 867         div1->n_col = div2->n_col;
 868
 869         equal = isl_mat_is_equal(div1, div2);
 870
 871         div1->n_row = n_row;
 872         div1->n_col = n_col;
 873
 874         return equal;
 875 }
 876
 877 static void expand_row(__isl_keep isl_mat *dst, int d,
 878         __isl_keep isl_mat *src, int s, int *exp)
 879 {
 880         int i;
 881         unsigned c = src->n_col - src->n_row;
 882
 883         isl_seq_cpy(dst->row[d], src->row[s], c);
 884         isl_seq_clr(dst->row[d] + c, dst->n_col - c);
 885
 886         for (i = 0; i < s; ++i)
 887                 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
 888 }
 889
 890 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
 891 {
 892         int li, lj;
 893
 894         li = isl_seq_last_non_zero(div->row[i], div->n_col);
 895         lj = isl_seq_last_non_zero(div->row[j], div->n_col);
 896
 897         if (li != lj)
 898                 return li - lj;
 899
 900         return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
 901 }
 902
 903 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
 904         __isl_keep isl_mat *div2, int *exp1, int *exp2)
 905 {
 906         int i, j, k;
 907         isl_mat *div = NULL;
 908         unsigned d = div1->n_col - div1->n_row;
 909
 910         div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
 911                                 d + div1->n_row + div2->n_row);
 912         if (!div)
 913                 return NULL;
 914
 915         for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
 916                 int cmp;
 917
 918                 expand_row(div, k, div1, i, exp1);
 919                 expand_row(div, k + 1, div2, j, exp2);
 920
 921                 cmp = cmp_row(div, k, k + 1);
 922                 if (cmp == 0) {
 923                         exp1[i++] = k;
 924                         exp2[j++] = k;
 925                 } else if (cmp < 0) {
 926                         exp1[i++] = k;
 927                 } else {
 928                         exp2[j++] = k;
 929                         isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
 930                 }
 931         }
 932         for (; i < div1->n_row; ++i, ++k) {
 933                 expand_row(div, k, div1, i, exp1);
 934                 exp1[i] = k;
 935         }
 936         for (; j < div2->n_row; ++j, ++k) {
 937                 expand_row(div, k, div2, j, exp2);
 938                 exp2[j] = k;
 939         }
 940
 941         div->n_row = k;
 942         div->n_col = d + k;
 943
 944         return div;
 945 }
 946
 947 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
 948         int *exp, int first)
 949 {
 950         int i;
 951         struct isl_upoly_rec *rec;
 952
 953         if (isl_upoly_is_cst(up))
 954                 return up;
 955
 956         if (up->var < first)
 957                 return up;
 958
 959         if (exp[up->var - first] == up->var - first)
 960                 return up;
 961
 962         up = isl_upoly_cow(up);
 963         if (!up)
 964                 goto error;
 965
 966         up->var = exp[up->var - first] + first;
 967
 968         rec = isl_upoly_as_rec(up);
 969         if (!rec)
 970                 goto error;
 971
 972         for (i = 0; i < rec->n; ++i) {
 973                 rec->p[i] = expand(rec->p[i], exp, first);
 974                 if (!rec->p[i])
 975                         goto error;
 976         }
 977
 978         return up;
 979 error:
 980         isl_upoly_free(up);
 981         return NULL;
 982 }
 983
 984 static __isl_give isl_qpolynomial *with_merged_divs(
 985         __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
 986                                           __isl_take isl_qpolynomial *qp2),
 987         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
 988 {
 989         int *exp1 = NULL;
 990         int *exp2 = NULL;
 991         isl_mat *div = NULL;
 992
 993         qp1 = isl_qpolynomial_cow(qp1);
 994         qp2 = isl_qpolynomial_cow(qp2);
 995
 996         if (!qp1 || !qp2)
 997                 goto error;
 998
 999         isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1000                                 qp1->div->n_col >= qp2->div->n_col, goto error);
1001
1002         exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1003         exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1004         if (!exp1 || !exp2)
1005                 goto error;
1006
1007         div = merge_divs(qp1->div, qp2->div, exp1, exp2);
1008         if (!div)
1009                 goto error;
1010
1011         isl_mat_free(qp1->div);
1012         qp1->div = isl_mat_copy(div);
1013         isl_mat_free(qp2->div);
1014         qp2->div = isl_mat_copy(div);
1015
1016         qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1017         qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1018
1019         if (!qp1->upoly || !qp2->upoly)
1020                 goto error;
1021
1022         isl_mat_free(div);
1023         free(exp1);
1024         free(exp2);
1025
1026         return fn(qp1, qp2);
1027 error:
1028         isl_mat_free(div);
1029         free(exp1);
1030         free(exp2);
1031         isl_qpolynomial_free(qp1);
1032         isl_qpolynomial_free(qp2);
1033         return NULL;
1034 }
1035
1036 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1037         __isl_take isl_qpolynomial *qp2)
1038 {
1039         qp1 = isl_qpolynomial_cow(qp1);
1040
1041         if (!qp1 || !qp2)
1042                 goto error;
1043
1044         if (qp1->div->n_row < qp2->div->n_row)
1045                 return isl_qpolynomial_add(qp2, qp1);
1046
1047         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1048         if (!compatible_divs(qp1->div, qp2->div))
1049                 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1050
1051         qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1052         if (!qp1->upoly)
1053                 goto error;
1054
1055         isl_qpolynomial_free(qp2);
1056
1057         return qp1;
1058 error:
1059         isl_qpolynomial_free(qp1);
1060         isl_qpolynomial_free(qp2);
1061         return NULL;
1062 }
1063
1064 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1065         __isl_take isl_qpolynomial *qp2)
1066 {
1067         return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1068 }
1069
1070 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1071 {
1072         qp = isl_qpolynomial_cow(qp);
1073
1074         if (!qp)
1075                 return NULL;
1076
1077         qp->upoly = isl_upoly_neg(qp->upoly);
1078         if (!qp->upoly)
1079                 goto error;
1080
1081         return qp;
1082 error:
1083         isl_qpolynomial_free(qp);
1084         return NULL;
1085 }
1086
1087 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1088         __isl_take isl_qpolynomial *qp2)
1089 {
1090         qp1 = isl_qpolynomial_cow(qp1);
1091
1092         if (!qp1 || !qp2)
1093                 goto error;
1094
1095         if (qp1->div->n_row < qp2->div->n_row)
1096                 return isl_qpolynomial_mul(qp2, qp1);
1097
1098         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1099         if (!compatible_divs(qp1->div, qp2->div))
1100                 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1101
1102         qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1103         if (!qp1->upoly)
1104                 goto error;
1105
1106         isl_qpolynomial_free(qp2);
1107
1108         return qp1;
1109 error:
1110         isl_qpolynomial_free(qp1);
1111         isl_qpolynomial_free(qp2);
1112         return NULL;
1113 }
1114
1115 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1116 {
1117         struct isl_qpolynomial *qp;
1118         struct isl_upoly_cst *cst;
1119
1120         qp = isl_qpolynomial_alloc(dim, 0);
1121         if (!qp)
1122                 return NULL;
1123
1124         qp->upoly = isl_upoly_zero(dim->ctx);
1125         if (!qp->upoly)
1126                 goto error;
1127
1128         return qp;
1129 error:
1130         isl_qpolynomial_free(qp);
1131         return NULL;
1132 }
1133
1134 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1135 {
1136         struct isl_qpolynomial *qp;
1137         struct isl_upoly_cst *cst;
1138
1139         qp = isl_qpolynomial_alloc(dim, 0);
1140         if (!qp)
1141                 return NULL;
1142
1143         qp->upoly = isl_upoly_infty(dim->ctx);
1144         if (!qp->upoly)
1145                 goto error;
1146
1147         return qp;
1148 error:
1149         isl_qpolynomial_free(qp);
1150         return NULL;
1151 }
1152
1153 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1154 {
1155         struct isl_qpolynomial *qp;
1156         struct isl_upoly_cst *cst;
1157
1158         qp = isl_qpolynomial_alloc(dim, 0);
1159         if (!qp)
1160                 return NULL;
1161
1162         qp->upoly = isl_upoly_nan(dim->ctx);
1163         if (!qp->upoly)
1164                 goto error;
1165
1166         return qp;
1167 error:
1168         isl_qpolynomial_free(qp);
1169         return NULL;
1170 }
1171
1172 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1173         isl_int v)
1174 {
1175         struct isl_qpolynomial *qp;
1176         struct isl_upoly_cst *cst;
1177
1178         qp = isl_qpolynomial_alloc(dim, 0);
1179         if (!qp)
1180                 return NULL;
1181
1182         qp->upoly = isl_upoly_zero(dim->ctx);
1183         if (!qp->upoly)
1184                 goto error;
1185         cst = isl_upoly_as_cst(qp->upoly);
1186         isl_int_set(cst->n, v);
1187
1188         return qp;
1189 error:
1190         isl_qpolynomial_free(qp);
1191         return NULL;
1192 }
1193
1194 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1195         isl_int *n, isl_int *d)
1196 {
1197         struct isl_upoly_cst *cst;
1198
1199         if (!qp)
1200                 return -1;
1201
1202         if (!isl_upoly_is_cst(qp->upoly))
1203                 return 0;
1204
1205         cst = isl_upoly_as_cst(qp->upoly);
1206         if (!cst)
1207                 return -1;
1208
1209         if (n)
1210                 isl_int_set(*n, cst->n);
1211         if (d)
1212                 isl_int_set(*d, cst->d);
1213
1214         return 1;
1215 }
1216
1217 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
1218         __isl_keep isl_qpolynomial *qp2)
1219 {
1220         if (!qp1 || !qp2)
1221                 return -1;
1222
1223         return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1224 }
1225
1226 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1227 {
1228         int i;
1229         struct isl_upoly_rec *rec;
1230
1231         if (isl_upoly_is_cst(up)) {
1232                 struct isl_upoly_cst *cst;
1233                 cst = isl_upoly_as_cst(up);
1234                 if (!cst)
1235                         return;
1236                 isl_int_lcm(*d, *d, cst->d);
1237                 return;
1238         }
1239
1240         rec = isl_upoly_as_rec(up);
1241         if (!rec)
1242                 return;
1243
1244         for (i = 0; i < rec->n; ++i)
1245                 upoly_update_den(rec->p[i], d);
1246 }
1247
1248 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1249 {
1250         isl_int_set_si(*d, 1);
1251         if (!qp)
1252                 return;
1253         upoly_update_den(qp->upoly, d);
1254 }
1255
1256 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1257         int pos, int power)
1258 {
1259         int i;
1260         struct isl_qpolynomial *qp;
1261         struct isl_upoly_rec *rec;
1262         struct isl_upoly_cst *cst;
1263         struct isl_ctx *ctx;
1264
1265         if (!dim)
1266                 return NULL;
1267
1268         ctx = dim->ctx;
1269
1270         qp = isl_qpolynomial_alloc(dim, 0);
1271         if (!qp)
1272                 return NULL;
1273
1274         qp->upoly = &isl_upoly_alloc_rec(ctx, pos, 1 + power)->up;
1275         if (!qp->upoly)
1276                 goto error;
1277         rec = isl_upoly_as_rec(qp->upoly);
1278         for (i = 0; i < 1 + power; ++i) {
1279                 rec->p[i] = isl_upoly_zero(ctx);
1280                 if (!rec->p[i])
1281                         goto error;
1282                 rec->n++;
1283         }
1284         cst = isl_upoly_as_cst(rec->p[power]);
1285         isl_int_set_si(cst->n, 1);
1286
1287         return qp;
1288 error:
1289         isl_qpolynomial_free(qp);
1290         return NULL;
1291 }
1292
1293 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1294         enum isl_dim_type type, unsigned pos)
1295 {
1296         if (!dim)
1297                 return NULL;
1298
1299         isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1300         isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1301
1302         if (type == isl_dim_set)
1303                 pos += isl_dim_size(dim, isl_dim_param);
1304
1305         return isl_qpolynomial_pow(dim, pos, 1);
1306 error:
1307         isl_dim_free(dim);
1308         return NULL;
1309 }
1310
1311 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1312         int power)
1313 {
1314         struct isl_qpolynomial *qp = NULL;
1315         struct isl_upoly_rec *rec;
1316         struct isl_upoly_cst *cst;
1317         int i;
1318         int pos;
1319
1320         if (!div)
1321                 return NULL;
1322         isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1323
1324         qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1);
1325         if (!qp)
1326                 goto error;
1327
1328         isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1329         isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1330
1331         pos = isl_dim_total(qp->dim);
1332         qp->upoly = &isl_upoly_alloc_rec(div->ctx, pos, 1 + power)->up;
1333         if (!qp->upoly)
1334                 goto error;
1335         rec = isl_upoly_as_rec(qp->upoly);
1336         for (i = 0; i < 1 + power; ++i) {
1337                 rec->p[i] = isl_upoly_zero(div->ctx);
1338                 if (!rec->p[i])
1339                         goto error;
1340                 rec->n++;
1341         }
1342         cst = isl_upoly_as_cst(rec->p[power]);
1343         isl_int_set_si(cst->n, 1);
1344
1345         isl_div_free(div);
1346
1347         return qp;
1348 error:
1349         isl_qpolynomial_free(qp);
1350         isl_div_free(div);
1351         return NULL;
1352 }
1353
1354 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1355 {
1356         return isl_qpolynomial_div_pow(div, 1);
1357 }
1358
1359 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1360         const isl_int n, const isl_int d)
1361 {
1362         struct isl_qpolynomial *qp;
1363         struct isl_upoly_cst *cst;
1364
1365         qp = isl_qpolynomial_alloc(dim, 0);
1366         if (!qp)
1367                 return NULL;
1368
1369         qp->upoly = isl_upoly_zero(dim->ctx);
1370         if (!qp->upoly)
1371                 goto error;
1372         cst = isl_upoly_as_cst(qp->upoly);
1373         isl_int_set(cst->n, n);
1374         isl_int_set(cst->d, d);
1375
1376         return qp;
1377 error:
1378         isl_qpolynomial_free(qp);
1379         return NULL;
1380 }
1381
1382 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
1383 {
1384         struct isl_upoly_rec *rec;
1385         int i;
1386
1387         if (!up)
1388                 return -1;
1389
1390         if (isl_upoly_is_cst(up))
1391                 return 0;
1392
1393         if (up->var < d)
1394                 active[up->var] = 1;
1395
1396         rec = isl_upoly_as_rec(up);
1397         for (i = 0; i < rec->n; ++i)
1398                 if (up_set_active(rec->p[i], active, d) < 0)
1399                         return -1;
1400
1401         return 0;
1402 }
1403
1404 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
1405 {
1406         int i, j;
1407         int d = isl_dim_total(qp->dim);
1408
1409         if (!qp || !active)
1410                 return -1;
1411
1412         for (i = 0; i < d; ++i)
1413                 for (j = 0; j < qp->div->n_row; ++j) {
1414                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
1415                                 continue;
1416                         active[i] = 1;
1417                         break;
1418                 }
1419
1420         return up_set_active(qp->upoly, active, d);
1421 }
1422
1423 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
1424         enum isl_dim_type type, unsigned first, unsigned n)
1425 {
1426         int i;
1427         int *active = NULL;
1428         int involves = 0;
1429
1430         if (!qp)
1431                 return -1;
1432         if (n == 0)
1433                 return 0;
1434
1435         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1436                         return -1);
1437         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1438                                  type == isl_dim_set, return -1);
1439
1440         active = isl_calloc_array(set->ctx, int, isl_dim_total(qp->dim));
1441         if (set_active(qp, active) < 0)
1442                 goto error;
1443
1444         if (type == isl_dim_set)
1445                 first += isl_dim_size(qp->dim, isl_dim_param);
1446         for (i = 0; i < n; ++i)
1447                 if (active[first + i]) {
1448                         involves = 1;
1449                         break;
1450                 }
1451
1452         free(active);
1453
1454         return involves;
1455 error:
1456         free(active);
1457         return -1;
1458 }
1459
1460 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
1461         enum isl_dim_type type, unsigned first, unsigned n)
1462 {
1463         int i;
1464
1465         if (!fold)
1466                 return -1;
1467         if (fold->n == 0 || n == 0)
1468                 return 0;
1469
1470         for (i = 0; i < fold->n; ++i) {
1471                 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
1472                                                             type, first, n);
1473                 if (involves < 0 || involves)
1474                         return involves;
1475         }
1476         return 0;
1477 }
1478
1479 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
1480         unsigned first, unsigned n)
1481 {
1482         int i;
1483         struct isl_upoly_rec *rec;
1484
1485         if (!up)
1486                 return NULL;
1487         if (n == 0 || up->var < 0 || up->var < first)
1488                 return up;
1489         if (up->var < first + n) {
1490                 up = replace_by_constant_term(up);
1491                 return isl_upoly_drop(up, first, n);
1492         }
1493         up = isl_upoly_cow(up);
1494         if (!up)
1495                 return NULL;
1496         up->var -= n;
1497         rec = isl_upoly_as_rec(up);
1498         if (!rec)
1499                 goto error;
1500
1501         for (i = 0; i < rec->n; ++i) {
1502                 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
1503                 if (!rec->p[i])
1504                         goto error;
1505         }
1506
1507         return up;
1508 error:
1509         isl_upoly_free(up);
1510         return NULL;
1511 }
1512
1513 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
1514         __isl_take isl_qpolynomial *qp,
1515         enum isl_dim_type type, unsigned first, unsigned n)
1516 {
1517         if (!qp)
1518                 return NULL;
1519         if (n == 0)
1520                 return qp;
1521
1522         qp = isl_qpolynomial_cow(qp);
1523         if (!qp)
1524                 return NULL;
1525
1526         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1527                         goto error);
1528         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1529                                  type == isl_dim_set, goto error);
1530
1531         qp->dim = isl_dim_drop(qp->dim, type, first, n);
1532         if (!qp->dim)
1533                 goto error;
1534
1535         if (type == isl_dim_set)
1536                 first += isl_dim_size(qp->dim, isl_dim_param);
1537
1538         qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
1539         if (!qp->div)
1540                 goto error;
1541
1542         qp->upoly = isl_upoly_drop(qp->upoly, first, n);
1543         if (!qp->upoly)
1544                 goto error;
1545
1546         return qp;
1547 error:
1548         isl_qpolynomial_free(qp);
1549         return NULL;
1550 }
1551
1552 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
1553         __isl_take isl_qpolynomial_fold *fold,
1554         enum isl_dim_type type, unsigned first, unsigned n)
1555 {
1556         int i;
1557
1558         if (!fold)
1559                 return NULL;
1560         if (n == 0)
1561                 return fold;
1562
1563         fold = isl_qpolynomial_fold_cow(fold);
1564         if (!fold)
1565                 return NULL;
1566         fold->dim = isl_dim_drop(fold->dim, type, first, n);
1567         if (!fold->dim)
1568                 goto error;
1569
1570         for (i = 0; i < fold->n; ++i) {
1571                 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
1572                                                             type, first, n);
1573                 if (!fold->qp[i])
1574                         goto error;
1575         }
1576
1577         return fold;
1578 error:
1579         isl_qpolynomial_fold_free(fold);
1580         return NULL;
1581 }
1582
1583 #undef PW
1584 #define PW isl_pw_qpolynomial
1585 #undef EL
1586 #define EL isl_qpolynomial
1587 #undef IS_ZERO
1588 #define IS_ZERO is_zero
1589 #undef FIELD
1590 #define FIELD qp
1591 #undef ADD
1592 #define ADD add
1593
1594 #include <isl_pw_templ.c>
1595
1596 #undef PW
1597 #define PW isl_pw_qpolynomial_fold
1598 #undef EL
1599 #define EL isl_qpolynomial_fold
1600 #undef IS_ZERO
1601 #define IS_ZERO is_empty
1602 #undef FIELD
1603 #define FIELD fold
1604 #undef ADD
1605 #define ADD fold
1606
1607 #include <isl_pw_templ.c>
1608
1609 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1610 {
1611         if (!pwqp)
1612                 return -1;
1613
1614         if (pwqp->n != -1)
1615                 return 0;
1616
1617         if (!isl_set_fast_is_universe(pwqp->p[0].set))
1618                 return 0;
1619
1620         return isl_qpolynomial_is_one(pwqp->p[0].qp);
1621 }
1622
1623 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1624         __isl_take isl_pw_qpolynomial *pwqp1,
1625         __isl_take isl_pw_qpolynomial *pwqp2)
1626 {
1627         int i, j, n;
1628         struct isl_pw_qpolynomial *res;
1629         isl_set *set;
1630
1631         if (!pwqp1 || !pwqp2)
1632                 goto error;
1633
1634         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1635                         goto error);
1636
1637         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1638                 isl_pw_qpolynomial_free(pwqp2);
1639                 return pwqp1;
1640         }
1641
1642         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1643                 isl_pw_qpolynomial_free(pwqp1);
1644                 return pwqp2;
1645         }
1646
1647         if (isl_pw_qpolynomial_is_one(pwqp1)) {
1648                 isl_pw_qpolynomial_free(pwqp1);
1649                 return pwqp2;
1650         }
1651
1652         if (isl_pw_qpolynomial_is_one(pwqp2)) {
1653                 isl_pw_qpolynomial_free(pwqp2);
1654                 return pwqp1;
1655         }
1656
1657         n = pwqp1->n * pwqp2->n;
1658         res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
1659
1660         for (i = 0; i < pwqp1->n; ++i) {
1661                 for (j = 0; j < pwqp2->n; ++j) {
1662                         struct isl_set *common;
1663                         struct isl_qpolynomial *prod;
1664                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1665                                                 isl_set_copy(pwqp2->p[j].set));
1666                         if (isl_set_fast_is_empty(common)) {
1667                                 isl_set_free(common);
1668                                 continue;
1669                         }
1670
1671                         prod = isl_qpolynomial_mul(
1672                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1673                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1674
1675                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
1676                 }
1677         }
1678
1679         isl_pw_qpolynomial_free(pwqp1);
1680         isl_pw_qpolynomial_free(pwqp2);
1681
1682         return res;
1683 error:
1684         isl_pw_qpolynomial_free(pwqp1);
1685         isl_pw_qpolynomial_free(pwqp2);
1686         return NULL;
1687 }
1688
1689 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1690         __isl_take isl_pw_qpolynomial *pwqp)
1691 {
1692         int i, j, n;
1693         struct isl_pw_qpolynomial *res;
1694         isl_set *set;
1695
1696         if (!pwqp)
1697                 return NULL;
1698
1699         if (isl_pw_qpolynomial_is_zero(pwqp))
1700                 return pwqp;
1701
1702         pwqp = isl_pw_qpolynomial_cow(pwqp);
1703
1704         for (i = 0; i < pwqp->n; ++i) {
1705                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1706                 if (!pwqp->p[i].qp)
1707                         goto error;
1708         }
1709
1710         return pwqp;
1711 error:
1712         isl_pw_qpolynomial_free(pwqp);
1713         return NULL;
1714 }
1715
1716 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1717         __isl_take isl_pw_qpolynomial *pwqp1,
1718         __isl_take isl_pw_qpolynomial *pwqp2)
1719 {
1720         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1721 }
1722
1723 __isl_give struct isl_upoly *isl_upoly_eval(
1724         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1725 {
1726         int i;
1727         struct isl_upoly_rec *rec;
1728         struct isl_upoly *res;
1729         struct isl_upoly *base;
1730
1731         if (isl_upoly_is_cst(up)) {
1732                 isl_vec_free(vec);
1733                 return up;
1734         }
1735
1736         rec = isl_upoly_as_rec(up);
1737         if (!rec)
1738                 goto error;
1739
1740         isl_assert(up->ctx, rec->n >= 1, goto error);
1741
1742         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1743
1744         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1745                                 isl_vec_copy(vec));
1746
1747         for (i = rec->n - 2; i >= 0; --i) {
1748                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1749                 res = isl_upoly_sum(res,
1750                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1751                                                             isl_vec_copy(vec)));
1752         }
1753
1754         isl_upoly_free(base);
1755         isl_upoly_free(up);
1756         isl_vec_free(vec);
1757         return res;
1758 error:
1759         isl_upoly_free(up);
1760         isl_vec_free(vec);
1761         return NULL;
1762 }
1763
1764 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1765         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1766 {
1767         isl_vec *ext;
1768         struct isl_upoly *up;
1769         isl_dim *dim;
1770
1771         if (!qp || !pnt)
1772                 goto error;
1773         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1774
1775         if (qp->div->n_row == 0)
1776                 ext = isl_vec_copy(pnt->vec);
1777         else {
1778                 int i;
1779                 unsigned dim = isl_dim_total(qp->dim);
1780                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1781                 if (!ext)
1782                         goto error;
1783
1784                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1785                 for (i = 0; i < qp->div->n_row; ++i) {
1786                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1787                                                 1 + dim + i, &ext->el[1+dim+i]);
1788                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1789                                         qp->div->row[i][0]);
1790                 }
1791         }
1792
1793         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1794         if (!up)
1795                 goto error;
1796
1797         dim = isl_dim_copy(qp->dim);
1798         isl_qpolynomial_free(qp);
1799         isl_point_free(pnt);
1800
1801         qp = isl_qpolynomial_alloc(dim, 0);
1802         if (!qp)
1803                 isl_upoly_free(up);
1804         else
1805                 qp->upoly = up;
1806
1807         return qp;
1808 error:
1809         isl_qpolynomial_free(qp);
1810         isl_point_free(pnt);
1811         return NULL;
1812 }
1813
1814 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
1815         __isl_keep struct isl_upoly_cst *cst2)
1816 {
1817         int cmp;
1818         isl_int t;
1819         isl_int_init(t);
1820         isl_int_mul(t, cst1->n, cst2->d);
1821         isl_int_submul(t, cst2->n, cst1->d);
1822         cmp = isl_int_sgn(t);
1823         isl_int_clear(t);
1824         return cmp;
1825 }
1826
1827 static __isl_give isl_qpolynomial *qpolynomial_min(
1828         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1829 {
1830         struct isl_upoly_cst *cst1, *cst2;
1831         int cmp;
1832
1833         if (!qp1 || !qp2)
1834                 goto error;
1835         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1836         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1837         cst1 = isl_upoly_as_cst(qp1->upoly);
1838         cst2 = isl_upoly_as_cst(qp2->upoly);
1839         cmp = isl_upoly_cmp(cst1, cst2);
1840
1841         if (cmp <= 0) {
1842                 isl_qpolynomial_free(qp2);
1843         } else {
1844                 isl_qpolynomial_free(qp1);
1845                 qp1 = qp2;
1846         }
1847         return qp1;
1848 error:
1849         isl_qpolynomial_free(qp1);
1850         isl_qpolynomial_free(qp2);
1851         return NULL;
1852 }
1853
1854 static __isl_give isl_qpolynomial *qpolynomial_max(
1855         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1856 {
1857         struct isl_upoly_cst *cst1, *cst2;
1858         int cmp;
1859
1860         if (!qp1 || !qp2)
1861                 goto error;
1862         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1863         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1864         cst1 = isl_upoly_as_cst(qp1->upoly);
1865         cst2 = isl_upoly_as_cst(qp2->upoly);
1866         cmp = isl_upoly_cmp(cst1, cst2);
1867
1868         if (cmp >= 0) {
1869                 isl_qpolynomial_free(qp2);
1870         } else {
1871                 isl_qpolynomial_free(qp1);
1872                 qp1 = qp2;
1873         }
1874         return qp1;
1875 error:
1876         isl_qpolynomial_free(qp1);
1877         isl_qpolynomial_free(qp2);
1878         return NULL;
1879 }
1880
1881 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
1882         __isl_keep isl_qpolynomial_fold *fold2)
1883 {
1884         int i;
1885
1886         if (!fold1 || !fold2)
1887                 return -1;
1888
1889         if (fold1->n != fold2->n)
1890                 return 0;
1891
1892         /* We probably want to sort the qps first... */
1893         for (i = 0; i < fold1->n; ++i) {
1894                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
1895                 if (eq < 0 || !eq)
1896                         return eq;
1897         }
1898
1899         return 1;
1900 }
1901
1902 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
1903         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
1904 {
1905         isl_qpolynomial *qp;
1906
1907         if (!fold || !pnt)
1908                 goto error;
1909         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
1910         isl_assert(pnt->dim->ctx,
1911                 fold->type == isl_fold_max || fold->type == isl_fold_min,
1912                 goto error);
1913
1914         if (fold->n == 0)
1915                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
1916         else {
1917                 int i;
1918                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
1919                                                 isl_point_copy(pnt));
1920                 for (i = 1; i < fold->n; ++i) {
1921                         isl_qpolynomial *qp_i;
1922                         qp_i = isl_qpolynomial_eval(
1923                                             isl_qpolynomial_copy(fold->qp[i]),
1924                                             isl_point_copy(pnt));
1925                         if (fold->type == isl_fold_max)
1926                                 qp = qpolynomial_max(qp, qp_i);
1927                         else
1928                                 qp = qpolynomial_min(qp, qp_i);
1929                 }
1930         }
1931         isl_qpolynomial_fold_free(fold);
1932         isl_point_free(pnt);
1933
1934         return qp;
1935 error:
1936         isl_qpolynomial_fold_free(fold);
1937         isl_point_free(pnt);
1938         return NULL;
1939 }
1940
1941 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1942         enum isl_dim_type dst_type, unsigned dst_pos,
1943         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1944 {
1945         if (!qp)
1946                 return NULL;
1947
1948         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1949         if (!qp->dim)
1950                 goto error;
1951
1952         /* Do something to polynomials when needed; later */
1953
1954         return qp;
1955 error:
1956         isl_qpolynomial_free(qp);
1957         return NULL;
1958 }
1959
1960 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1961         __isl_take isl_pw_qpolynomial *pwqp,
1962         enum isl_dim_type dst_type, unsigned dst_pos,
1963         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1964 {
1965         int i;
1966
1967         if (!pwqp)
1968                 return NULL;
1969
1970         pwqp->dim = isl_dim_move(pwqp->dim,
1971                                     dst_type, dst_pos, src_type, src_pos, n);
1972         if (!pwqp->dim)
1973                 goto error;
1974
1975         for (i = 0; i < pwqp->n; ++i) {
1976                 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1977                                                 src_type, src_pos, n);
1978                 if (!pwqp->p[i].set)
1979                         goto error;
1980                 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1981                                         dst_type, dst_pos, src_type, src_pos, n);
1982                 if (!pwqp->p[i].qp)
1983                         goto error;
1984         }
1985
1986         return pwqp;
1987 error:
1988         isl_pw_qpolynomial_free(pwqp);
1989         return NULL;
1990 }
1991
1992 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
1993         __isl_keep isl_pw_qpolynomial *pwqp)
1994 {
1995         if (!pwqp)
1996                 return NULL;
1997
1998         return isl_dim_copy(pwqp->dim);
1999 }
2000
2001 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
2002         enum isl_dim_type type)
2003 {
2004         return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
2005 }
2006
2007 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
2008         __isl_take isl_mat *div)
2009 {
2010         isl_term *term;
2011         int n;
2012
2013         if (!dim || !div)
2014                 goto error;
2015
2016         n = isl_dim_total(dim) + div->n_row;
2017
2018         term = isl_calloc(dim->ctx, struct isl_term,
2019                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
2020         if (!term)
2021                 goto error;
2022
2023         term->ref = 1;
2024         term->dim = dim;
2025         term->div = div;
2026         isl_int_init(term->n);
2027         isl_int_init(term->d);
2028
2029         return term;
2030 error:
2031         isl_dim_free(dim);
2032         isl_mat_free(div);
2033         return NULL;
2034 }
2035
2036 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
2037 {
2038         if (!term)
2039                 return NULL;
2040
2041         term->ref++;
2042         return term;
2043 }
2044
2045 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
2046 {
2047         int i;
2048         isl_term *dup;
2049         unsigned total;
2050
2051         if (term)
2052                 return NULL;
2053
2054         total = isl_dim_total(term->dim) + term->div->n_row;
2055
2056         dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
2057         if (!dup)
2058                 return NULL;
2059
2060         isl_int_set(dup->n, term->n);
2061         isl_int_set(dup->d, term->d);
2062
2063         for (i = 0; i < total; ++i)
2064                 dup->pow[i] = term->pow[i];
2065
2066         return dup;
2067 }
2068
2069 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
2070 {
2071         if (!term)
2072                 return NULL;
2073
2074         if (term->ref == 1)
2075                 return term;
2076         term->ref--;
2077         return isl_term_dup(term);
2078 }
2079
2080 void isl_term_free(__isl_take isl_term *term)
2081 {
2082         if (!term)
2083                 return;
2084
2085         if (--term->ref > 0)
2086                 return;
2087
2088         isl_dim_free(term->dim);
2089         isl_mat_free(term->div);
2090         isl_int_clear(term->n);
2091         isl_int_clear(term->d);
2092         free(term);
2093 }
2094
2095 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
2096 {
2097         if (!term)
2098                 return 0;
2099
2100         switch (type) {
2101         case isl_dim_param:
2102         case isl_dim_in:
2103         case isl_dim_out:       return isl_dim_size(term->dim, type);
2104         case isl_dim_div:       return term->div->n_row;
2105         case isl_dim_all:       return isl_dim_total(term->dim) + term->div->n_row;
2106         }
2107 }
2108
2109 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
2110 {
2111         return term ? term->dim->ctx : NULL;
2112 }
2113
2114 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
2115 {
2116         if (!term)
2117                 return;
2118         isl_int_set(*n, term->n);
2119 }
2120
2121 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
2122 {
2123         if (!term)
2124                 return;
2125         isl_int_set(*d, term->d);
2126 }
2127
2128 int isl_term_get_exp(__isl_keep isl_term *term,
2129         enum isl_dim_type type, unsigned pos)
2130 {
2131         if (!term)
2132                 return -1;
2133
2134         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
2135
2136         if (type >= isl_dim_set)
2137                 pos += isl_dim_size(term->dim, isl_dim_param);
2138         if (type >= isl_dim_div)
2139                 pos += isl_dim_size(term->dim, isl_dim_set);
2140
2141         return term->pow[pos];
2142 }
2143
2144 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
2145 {
2146         isl_basic_map *bmap;
2147         unsigned total;
2148         int k;
2149
2150         if (!term)
2151                 return NULL;
2152
2153         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
2154                         return NULL);
2155
2156         total = term->div->n_col - term->div->n_row - 2;
2157         /* No nested divs for now */
2158         isl_assert(term->dim->ctx,
2159                 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
2160                                         term->div->n_row) == -1,
2161                 return NULL);
2162
2163         bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
2164         if ((k = isl_basic_map_alloc_div(bmap)) < 0)
2165                 goto error;
2166
2167         isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
2168
2169         return isl_basic_map_div(bmap, k);
2170 error:
2171         isl_basic_map_free(bmap);
2172         return NULL;
2173 }
2174
2175 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
2176         int (*fn)(__isl_take isl_term *term, void *user),
2177         __isl_take isl_term *term, void *user)
2178 {
2179         int i;
2180         struct isl_upoly_rec *rec;
2181
2182         if (!up || !term)
2183                 goto error;
2184
2185         if (isl_upoly_is_zero(up))
2186                 return term;
2187
2188         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
2189         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
2190         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
2191
2192         if (isl_upoly_is_cst(up)) {
2193                 struct isl_upoly_cst *cst;
2194                 cst = isl_upoly_as_cst(up);
2195                 if (!cst)
2196                         goto error;
2197                 term = isl_term_cow(term);
2198                 if (!term)
2199                         goto error;
2200                 isl_int_set(term->n, cst->n);
2201                 isl_int_set(term->d, cst->d);
2202                 if (fn(isl_term_copy(term), user) < 0)
2203                         goto error;
2204                 return term;
2205         }
2206
2207         rec = isl_upoly_as_rec(up);
2208         if (!rec)
2209                 goto error;
2210
2211         for (i = 0; i < rec->n; ++i) {
2212                 term = isl_term_cow(term);
2213                 if (!term)
2214                         goto error;
2215                 term->pow[up->var] = i;
2216                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2217                 if (!term)
2218                         goto error;
2219         }
2220         term->pow[up->var] = 0;
2221
2222         return term;
2223 error:
2224         isl_term_free(term);
2225         return NULL;
2226 }
2227
2228 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
2229         int (*fn)(__isl_take isl_term *term, void *user), void *user)
2230 {
2231         isl_term *term;
2232
2233         if (!qp)
2234                 return -1;
2235
2236         term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
2237         if (!term)
2238                 return -1;
2239
2240         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
2241
2242         isl_term_free(term);
2243
2244         return term ? 0 : -1;
2245 }
2246
2247 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2248         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2249                     void *user), void *user)
2250 {
2251         int i;
2252
2253         if (!pwqp)
2254                 return -1;
2255
2256         for (i = 0; i < pwqp->n; ++i)
2257                 if (fn(isl_set_copy(pwqp->p[i].set),
2258                                 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
2259                         return -1;
2260
2261         return 0;
2262 }
2263
2264 static int any_divs(__isl_keep isl_set *set)
2265 {
2266         int i;
2267
2268         if (!set)
2269                 return -1;
2270
2271         for (i = 0; i < set->n; ++i)
2272                 if (set->p[i]->n_div > 0)
2273                         return 1;
2274
2275         return 0;
2276 }
2277
2278 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
2279         __isl_take isl_dim *dim)
2280 {
2281         if (!qp || !dim)
2282                 goto error;
2283
2284         if (isl_dim_equal(qp->dim, dim)) {
2285                 isl_dim_free(dim);
2286                 return qp;
2287         }
2288
2289         qp = isl_qpolynomial_cow(qp);
2290         if (!qp)
2291                 goto error;
2292
2293         if (qp->div->n_row) {
2294                 int i;
2295                 int extra;
2296                 unsigned total;
2297                 int *exp;
2298
2299                 extra = isl_dim_size(dim, isl_dim_set) -
2300                                 isl_dim_size(qp->dim, isl_dim_set);
2301                 total = isl_dim_total(qp->dim);
2302                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
2303                 if (!exp)
2304                         goto error;
2305                 for (i = 0; i < qp->div->n_row; ++i)
2306                         exp[i] = extra + i;
2307                 qp->upoly = expand(qp->upoly, exp, total);
2308                 free(exp);
2309                 if (!qp->upoly)
2310                         goto error;
2311                 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
2312                 if (!qp->div)
2313                         goto error;
2314                 for (i = 0; i < qp->div->n_row; ++i)
2315                         isl_seq_clr(qp->div->row[i] + 2 + total, extra);
2316         }
2317
2318         isl_dim_free(qp->dim);
2319         qp->dim = dim;
2320
2321         return qp;
2322 error:
2323         isl_dim_free(dim);
2324         isl_qpolynomial_free(qp);
2325         return NULL;
2326 }
2327
2328 static int foreach_lifted_subset(__isl_take isl_set *set,
2329         __isl_take isl_qpolynomial *qp,
2330         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2331                     void *user), void *user)
2332 {
2333         int i;
2334
2335         if (!set || !qp)
2336                 goto error;
2337
2338         for (i = 0; i < set->n; ++i) {
2339                 isl_set *lift;
2340                 isl_qpolynomial *copy;
2341
2342                 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2343                 lift = isl_set_lift(lift);
2344
2345                 copy = isl_qpolynomial_copy(qp);
2346                 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2347
2348                 if (fn(lift, copy, user) < 0)
2349                         goto error;
2350         }
2351
2352         isl_set_free(set);
2353         isl_qpolynomial_free(qp);
2354
2355         return 0;
2356 error:
2357         isl_set_free(set);
2358         isl_qpolynomial_free(qp);
2359         return -1;
2360 }
2361
2362 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2363         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2364                     void *user), void *user)
2365 {
2366         int i;
2367
2368         if (!pwqp)
2369                 return -1;
2370
2371         for (i = 0; i < pwqp->n; ++i) {
2372                 isl_set *set;
2373                 isl_qpolynomial *qp;
2374
2375                 set = isl_set_copy(pwqp->p[i].set);
2376                 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2377                 if (!any_divs(set)) {
2378                         if (fn(set, qp, user) < 0)
2379                                 return -1;
2380                         continue;
2381                 }
2382                 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2383                         return -1;
2384         }
2385
2386         return 0;
2387 }
2388
2389 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
2390         enum isl_fold type, __isl_take isl_dim *dim, int n)
2391 {
2392         isl_qpolynomial_fold *fold;
2393
2394         if (!dim)
2395                 goto error;
2396
2397         isl_assert(dim->ctx, n >= 0, goto error);
2398         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
2399                         sizeof(struct isl_qpolynomial_fold) +
2400                         (n - 1) * sizeof(struct isl_qpolynomial *));
2401         if (!fold)
2402                 goto error;
2403
2404         fold->ref = 1;
2405         fold->size = n;
2406         fold->n = 0;
2407         fold->type = type;
2408         fold->dim = dim;
2409
2410         return fold;
2411 error:
2412         isl_dim_free(dim);
2413         return NULL;
2414 }
2415
2416 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
2417         __isl_take isl_dim *dim)
2418 {
2419         return qpolynomial_fold_alloc(type, dim, 0);
2420 }
2421
2422 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
2423         enum isl_fold type, __isl_take isl_qpolynomial *qp)
2424 {
2425         isl_qpolynomial_fold *fold;
2426
2427         if (!qp)
2428                 return NULL;
2429
2430         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
2431         if (!fold)
2432                 goto error;
2433
2434         fold->qp[0] = qp;
2435         fold->n++;
2436
2437         return fold;
2438 error:
2439         isl_qpolynomial_fold_free(fold);
2440         isl_qpolynomial_free(qp);
2441         return NULL;
2442 }
2443
2444 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
2445         __isl_keep isl_qpolynomial_fold *fold)
2446 {
2447         if (!fold)
2448                 return NULL;
2449
2450         fold->ref++;
2451         return fold;
2452 }
2453
2454 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
2455         __isl_keep isl_qpolynomial_fold *fold)
2456 {
2457         int i;
2458         isl_qpolynomial_fold *dup;
2459
2460         if (!fold)
2461                 return NULL;
2462         dup = qpolynomial_fold_alloc(fold->type,
2463                                         isl_dim_copy(fold->dim), fold->n);
2464         if (!dup)
2465                 return NULL;
2466
2467         for (i = 0; i < fold->n; ++i) {
2468                 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
2469                 if (!dup->qp[i])
2470                         goto error;
2471         }
2472
2473         return dup;
2474 error:
2475         isl_qpolynomial_fold_free(dup);
2476         return NULL;
2477 }
2478
2479 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
2480         __isl_take isl_qpolynomial_fold *fold)
2481 {
2482         if (!fold)
2483                 return NULL;
2484
2485         if (fold->ref == 1)
2486                 return fold;
2487         fold->ref--;
2488         return isl_qpolynomial_fold_dup(fold);
2489 }
2490
2491 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
2492 {
2493         int i;
2494
2495         if (!fold)
2496                 return;
2497         if (--fold->ref > 0)
2498                 return;
2499
2500         for (i = 0; i < fold->n; ++i)
2501                 isl_qpolynomial_free(fold->qp[i]);
2502         isl_dim_free(fold->dim);
2503         free(fold);
2504 }
2505
2506 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
2507 {
2508         if (!fold)
2509                 return -1;
2510
2511         return fold->n == 0;
2512 }
2513
2514 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
2515         __isl_take isl_qpolynomial_fold *fold1,
2516         __isl_take isl_qpolynomial_fold *fold2)
2517 {
2518         int i;
2519         struct isl_qpolynomial_fold *res = NULL;
2520
2521         if (!fold1 || !fold2)
2522                 goto error;
2523
2524         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
2525         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
2526                         goto error);
2527
2528         if (isl_qpolynomial_fold_is_empty(fold1)) {
2529                 isl_qpolynomial_fold_free(fold1);
2530                 return fold2;
2531         }
2532
2533         if (isl_qpolynomial_fold_is_empty(fold2)) {
2534                 isl_qpolynomial_fold_free(fold2);
2535                 return fold1;
2536         }
2537
2538         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
2539                                         fold1->n + fold2->n);
2540         if (!res)
2541                 goto error;
2542
2543         for (i = 0; i < fold1->n; ++i) {
2544                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
2545                 if (!res->qp[res->n])
2546                         goto error;
2547                 res->n++;
2548         }
2549
2550         for (i = 0; i < fold2->n; ++i) {
2551                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
2552                 if (!res->qp[res->n])
2553                         goto error;
2554                 res->n++;
2555         }
2556
2557         isl_qpolynomial_fold_free(fold1);
2558         isl_qpolynomial_fold_free(fold2);
2559
2560         return res;
2561 error:
2562         isl_qpolynomial_fold_free(res);
2563         isl_qpolynomial_fold_free(fold1);
2564         isl_qpolynomial_fold_free(fold2);
2565         return NULL;
2566 }
2567
2568 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2569         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
2570 {
2571         int i;
2572         isl_pw_qpolynomial_fold *pwf;
2573
2574         if (!pwqp)
2575                 return NULL;
2576
2577         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
2578
2579         for (i = 0; i < pwqp->n; ++i)
2580                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
2581                         isl_set_copy(pwqp->p[i].set),
2582                         isl_qpolynomial_fold_alloc(type,
2583                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
2584
2585         isl_pw_qpolynomial_free(pwqp);
2586
2587         return pwf;
2588 }
2589
2590 /* For each parameter or variable that does not appear in qp,
2591  * first eliminate the variable from all constraints and then set it to zero.
2592  */
2593 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
2594         __isl_keep isl_qpolynomial *qp)
2595 {
2596         int *active = NULL;
2597         int i;
2598         int d;
2599         unsigned nparam;
2600         unsigned nvar;
2601
2602         if (!set || !qp)
2603                 goto error;
2604
2605         d = isl_dim_total(set->dim);
2606         active = isl_calloc_array(set->ctx, int, d);
2607         if (set_active(qp, active) < 0)
2608                 goto error;
2609
2610         for (i = 0; i < d; ++i)
2611                 if (!active[i])
2612                         break;
2613
2614         if (i == d) {
2615                 free(active);
2616                 return set;
2617         }
2618
2619         nparam = isl_dim_size(set->dim, isl_dim_param);
2620         nvar = isl_dim_size(set->dim, isl_dim_set);
2621         for (i = 0; i < nparam; ++i) {
2622                 if (active[i])
2623                         continue;
2624                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
2625                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
2626         }
2627         for (i = 0; i < nvar; ++i) {
2628                 if (active[i])
2629                         continue;
2630                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
2631                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
2632         }
2633
2634         free(active);
2635
2636         return set;
2637 error:
2638         free(active);
2639         isl_set_free(set);
2640         return NULL;
2641 }
2642
2643 struct isl_max_data {
2644         isl_qpolynomial *qp;
2645         int first;
2646         isl_qpolynomial *max;
2647 };
2648
2649 static int max_fn(__isl_take isl_point *pnt, void *user)
2650 {
2651         struct isl_max_data *data = (struct isl_max_data *)user;
2652         isl_qpolynomial *val;
2653
2654         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
2655         if (data->first) {
2656                 data->first = 0;
2657                 data->max = val;
2658         } else {
2659                 data->max = qpolynomial_max(data->max, val);
2660         }
2661
2662         return 0;
2663 }
2664
2665 static __isl_give isl_qpolynomial *guarded_qpolynomial_max(
2666         __isl_take isl_set *set, __isl_take isl_qpolynomial *qp)
2667 {
2668         struct isl_max_data data = { NULL, 1, NULL };
2669
2670         if (!set || !qp)
2671                 goto error;
2672
2673         if (isl_upoly_is_cst(qp->upoly)) {
2674                 isl_set_free(set);
2675                 return qp;
2676         }
2677
2678         set = fix_inactive(set, qp);
2679
2680         data.qp = qp;
2681         if (isl_set_foreach_point(set, max_fn, &data) < 0)
2682                 goto error;
2683
2684         isl_set_free(set);
2685         isl_qpolynomial_free(qp);
2686         return data.max;
2687 error:
2688         isl_set_free(set);
2689         isl_qpolynomial_free(qp);
2690         isl_qpolynomial_free(data.max);
2691         return NULL;
2692 }
2693
2694 /* Compute the maximal value attained by the piecewise quasipolynomial
2695  * on its domain or zero if the domain is empty.
2696  * In the worst case, the domain is scanned completely,
2697  * so the domain is assumed to be bounded.
2698  */
2699 __isl_give isl_qpolynomial *isl_pw_qpolynomial_max(
2700         __isl_take isl_pw_qpolynomial *pwqp)
2701 {
2702         int i;
2703         isl_qpolynomial *max;
2704
2705         if (!pwqp)
2706                 return NULL;
2707
2708         if (pwqp->n == 0) {
2709                 isl_dim *dim = isl_dim_copy(pwqp->dim);
2710                 isl_pw_qpolynomial_free(pwqp);
2711                 return isl_qpolynomial_zero(dim);
2712         }
2713
2714         max = guarded_qpolynomial_max(isl_set_copy(pwqp->p[0].set),
2715                                         isl_qpolynomial_copy(pwqp->p[0].qp));
2716         for (i = 1; i < pwqp->n; ++i) {
2717                 isl_qpolynomial *max_i;
2718                 max_i = guarded_qpolynomial_max(isl_set_copy(pwqp->p[i].set),
2719                                             isl_qpolynomial_copy(pwqp->p[i].qp));
2720                 max = qpolynomial_max(max, max_i);
2721         }
2722
2723         isl_pw_qpolynomial_free(pwqp);
2724         return max;
2725 }