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