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_mul(
1323         __isl_take isl_pw_qpolynomial *pwqp1,
1324         __isl_take isl_pw_qpolynomial *pwqp2)
1325 {
1326         int i, j, n;
1327         struct isl_pw_qpolynomial *res;
1328         isl_set *set;
1329
1330         if (!pwqp1 || !pwqp2)
1331                 goto error;
1332
1333         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1334                         goto error);
1335
1336         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1337                 isl_pw_qpolynomial_free(pwqp2);
1338                 return pwqp1;
1339         }
1340
1341         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1342                 isl_pw_qpolynomial_free(pwqp1);
1343                 return pwqp2;
1344         }
1345
1346         if (isl_pw_qpolynomial_is_one(pwqp1)) {
1347                 isl_pw_qpolynomial_free(pwqp1);
1348                 return pwqp2;
1349         }
1350
1351         if (isl_pw_qpolynomial_is_one(pwqp2)) {
1352                 isl_pw_qpolynomial_free(pwqp2);
1353                 return pwqp1;
1354         }
1355
1356         n = pwqp1->n * pwqp2->n;
1357         res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
1358
1359         for (i = 0; i < pwqp1->n; ++i) {
1360                 for (j = 0; j < pwqp2->n; ++j) {
1361                         struct isl_set *common;
1362                         struct isl_qpolynomial *prod;
1363                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1364                                                 isl_set_copy(pwqp2->p[j].set));
1365                         if (isl_set_fast_is_empty(common)) {
1366                                 isl_set_free(common);
1367                                 continue;
1368                         }
1369
1370                         prod = isl_qpolynomial_mul(
1371                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1372                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1373
1374                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
1375                 }
1376         }
1377
1378         isl_pw_qpolynomial_free(pwqp1);
1379         isl_pw_qpolynomial_free(pwqp2);
1380
1381         return res;
1382 error:
1383         isl_pw_qpolynomial_free(pwqp1);
1384         isl_pw_qpolynomial_free(pwqp2);
1385         return NULL;
1386 }
1387
1388 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1389         __isl_take isl_pw_qpolynomial *pwqp)
1390 {
1391         int i, j, n;
1392         struct isl_pw_qpolynomial *res;
1393         isl_set *set;
1394
1395         if (!pwqp)
1396                 return NULL;
1397
1398         if (isl_pw_qpolynomial_is_zero(pwqp))
1399                 return pwqp;
1400
1401         pwqp = isl_pw_qpolynomial_cow(pwqp);
1402
1403         for (i = 0; i < pwqp->n; ++i) {
1404                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1405                 if (!pwqp->p[i].qp)
1406                         goto error;
1407         }
1408
1409         return pwqp;
1410 error:
1411         isl_pw_qpolynomial_free(pwqp);
1412         return NULL;
1413 }
1414
1415 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1416         __isl_take isl_pw_qpolynomial *pwqp1,
1417         __isl_take isl_pw_qpolynomial *pwqp2)
1418 {
1419         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1420 }
1421
1422 __isl_give struct isl_upoly *isl_upoly_eval(
1423         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1424 {
1425         int i;
1426         struct isl_upoly_rec *rec;
1427         struct isl_upoly *res;
1428         struct isl_upoly *base;
1429
1430         if (isl_upoly_is_cst(up)) {
1431                 isl_vec_free(vec);
1432                 return up;
1433         }
1434
1435         rec = isl_upoly_as_rec(up);
1436         if (!rec)
1437                 goto error;
1438
1439         isl_assert(up->ctx, rec->n >= 1, goto error);
1440
1441         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1442
1443         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1444                                 isl_vec_copy(vec));
1445
1446         for (i = rec->n - 2; i >= 0; --i) {
1447                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1448                 res = isl_upoly_sum(res,
1449                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1450                                                             isl_vec_copy(vec)));
1451         }
1452
1453         isl_upoly_free(base);
1454         isl_upoly_free(up);
1455         isl_vec_free(vec);
1456         return res;
1457 error:
1458         isl_upoly_free(up);
1459         isl_vec_free(vec);
1460         return NULL;
1461 }
1462
1463 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1464         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1465 {
1466         isl_vec *ext;
1467         struct isl_upoly *up;
1468         isl_dim *dim;
1469
1470         if (!qp || !pnt)
1471                 goto error;
1472         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1473
1474         if (qp->div->n_row == 0)
1475                 ext = isl_vec_copy(pnt->vec);
1476         else {
1477                 int i;
1478                 unsigned dim = isl_dim_total(qp->dim);
1479                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1480                 if (!ext)
1481                         goto error;
1482
1483                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1484                 for (i = 0; i < qp->div->n_row; ++i) {
1485                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1486                                                 1 + dim + i, &ext->el[1+dim+i]);
1487                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1488                                         qp->div->row[i][0]);
1489                 }
1490         }
1491
1492         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1493         if (!up)
1494                 goto error;
1495
1496         dim = isl_dim_copy(qp->dim);
1497         isl_qpolynomial_free(qp);
1498         isl_point_free(pnt);
1499
1500         qp = isl_qpolynomial_alloc(dim, 0);
1501         if (!qp)
1502                 isl_upoly_free(up);
1503         else
1504                 qp->upoly = up;
1505
1506         return qp;
1507 error:
1508         isl_qpolynomial_free(qp);
1509         isl_point_free(pnt);
1510         return NULL;
1511 }
1512
1513 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
1514         __isl_keep struct isl_upoly_cst *cst2)
1515 {
1516         int cmp;
1517         isl_int t;
1518         isl_int_init(t);
1519         isl_int_mul(t, cst1->n, cst2->d);
1520         isl_int_submul(t, cst2->n, cst1->d);
1521         cmp = isl_int_sgn(t);
1522         isl_int_clear(t);
1523         return cmp;
1524 }
1525
1526 static __isl_give isl_qpolynomial *qpolynomial_min(
1527         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1528 {
1529         struct isl_upoly_cst *cst1, *cst2;
1530         int cmp;
1531
1532         if (!qp1 || !qp2)
1533                 goto error;
1534         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1535         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1536         cst1 = isl_upoly_as_cst(qp1->upoly);
1537         cst2 = isl_upoly_as_cst(qp2->upoly);
1538         cmp = isl_upoly_cmp(cst1, cst2);
1539
1540         if (cmp <= 0) {
1541                 isl_qpolynomial_free(qp2);
1542         } else {
1543                 isl_qpolynomial_free(qp1);
1544                 qp1 = qp2;
1545         }
1546         return qp1;
1547 error:
1548         isl_qpolynomial_free(qp1);
1549         isl_qpolynomial_free(qp2);
1550         return NULL;
1551 }
1552
1553 static __isl_give isl_qpolynomial *qpolynomial_max(
1554         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1555 {
1556         struct isl_upoly_cst *cst1, *cst2;
1557         int cmp;
1558
1559         if (!qp1 || !qp2)
1560                 goto error;
1561         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1562         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1563         cst1 = isl_upoly_as_cst(qp1->upoly);
1564         cst2 = isl_upoly_as_cst(qp2->upoly);
1565         cmp = isl_upoly_cmp(cst1, cst2);
1566
1567         if (cmp >= 0) {
1568                 isl_qpolynomial_free(qp2);
1569         } else {
1570                 isl_qpolynomial_free(qp1);
1571                 qp1 = qp2;
1572         }
1573         return qp1;
1574 error:
1575         isl_qpolynomial_free(qp1);
1576         isl_qpolynomial_free(qp2);
1577         return NULL;
1578 }
1579
1580 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
1581         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
1582 {
1583         isl_qpolynomial *qp;
1584
1585         if (!fold || !pnt)
1586                 goto error;
1587         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
1588         isl_assert(pnt->dim->ctx,
1589                 fold->type == isl_fold_max || fold->type == isl_fold_min,
1590                 goto error);
1591
1592         if (fold->n == 0)
1593                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
1594         else {
1595                 int i;
1596                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
1597                                                 isl_point_copy(pnt));
1598                 for (i = 1; i < fold->n; ++i) {
1599                         isl_qpolynomial *qp_i;
1600                         qp_i = isl_qpolynomial_eval(
1601                                             isl_qpolynomial_copy(fold->qp[i]),
1602                                             isl_point_copy(pnt));
1603                         if (fold->type == isl_fold_max)
1604                                 qp = qpolynomial_max(qp, qp_i);
1605                         else
1606                                 qp = qpolynomial_min(qp, qp_i);
1607                 }
1608         }
1609         isl_qpolynomial_fold_free(fold);
1610         isl_point_free(pnt);
1611
1612         return qp;
1613 error:
1614         isl_qpolynomial_fold_free(fold);
1615         isl_point_free(pnt);
1616         return NULL;
1617 }
1618
1619 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1620         enum isl_dim_type dst_type, unsigned dst_pos,
1621         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1622 {
1623         if (!qp)
1624                 return NULL;
1625
1626         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1627         if (!qp->dim)
1628                 goto error;
1629
1630         /* Do something to polynomials when needed; later */
1631
1632         return qp;
1633 error:
1634         isl_qpolynomial_free(qp);
1635         return NULL;
1636 }
1637
1638 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1639         __isl_take isl_pw_qpolynomial *pwqp,
1640         enum isl_dim_type dst_type, unsigned dst_pos,
1641         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1642 {
1643         int i;
1644
1645         if (!pwqp)
1646                 return NULL;
1647
1648         pwqp->dim = isl_dim_move(pwqp->dim,
1649                                     dst_type, dst_pos, src_type, src_pos, n);
1650         if (!pwqp->dim)
1651                 goto error;
1652
1653         for (i = 0; i < pwqp->n; ++i) {
1654                 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1655                                                 src_type, src_pos, n);
1656                 if (!pwqp->p[i].set)
1657                         goto error;
1658                 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1659                                         dst_type, dst_pos, src_type, src_pos, n);
1660                 if (!pwqp->p[i].qp)
1661                         goto error;
1662         }
1663
1664         return pwqp;
1665 error:
1666         isl_pw_qpolynomial_free(pwqp);
1667         return NULL;
1668 }
1669
1670 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
1671         __isl_keep isl_pw_qpolynomial *pwqp)
1672 {
1673         if (!pwqp)
1674                 return NULL;
1675
1676         return isl_dim_copy(pwqp->dim);
1677 }
1678
1679 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
1680         enum isl_dim_type type)
1681 {
1682         return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
1683 }
1684
1685 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
1686         __isl_take isl_mat *div)
1687 {
1688         isl_term *term;
1689         int n;
1690
1691         if (!dim || !div)
1692                 goto error;
1693
1694         n = isl_dim_total(dim) + div->n_row;
1695
1696         term = isl_calloc(dim->ctx, struct isl_term,
1697                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
1698         if (!term)
1699                 goto error;
1700
1701         term->ref = 1;
1702         term->dim = dim;
1703         term->div = div;
1704         isl_int_init(term->n);
1705         isl_int_init(term->d);
1706
1707         return term;
1708 error:
1709         isl_dim_free(dim);
1710         isl_mat_free(div);
1711         return NULL;
1712 }
1713
1714 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
1715 {
1716         if (!term)
1717                 return NULL;
1718
1719         term->ref++;
1720         return term;
1721 }
1722
1723 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
1724 {
1725         int i;
1726         isl_term *dup;
1727         unsigned total;
1728
1729         if (term)
1730                 return NULL;
1731
1732         total = isl_dim_total(term->dim) + term->div->n_row;
1733
1734         dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
1735         if (!dup)
1736                 return NULL;
1737
1738         isl_int_set(dup->n, term->n);
1739         isl_int_set(dup->d, term->d);
1740
1741         for (i = 0; i < total; ++i)
1742                 dup->pow[i] = term->pow[i];
1743
1744         return dup;
1745 }
1746
1747 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
1748 {
1749         if (!term)
1750                 return NULL;
1751
1752         if (term->ref == 1)
1753                 return term;
1754         term->ref--;
1755         return isl_term_dup(term);
1756 }
1757
1758 void isl_term_free(__isl_take isl_term *term)
1759 {
1760         if (!term)
1761                 return;
1762
1763         if (--term->ref > 0)
1764                 return;
1765
1766         isl_dim_free(term->dim);
1767         isl_mat_free(term->div);
1768         isl_int_clear(term->n);
1769         isl_int_clear(term->d);
1770         free(term);
1771 }
1772
1773 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
1774 {
1775         if (!term)
1776                 return 0;
1777
1778         switch (type) {
1779         case isl_dim_param:
1780         case isl_dim_in:
1781         case isl_dim_out:       return isl_dim_size(term->dim, type);
1782         case isl_dim_div:       return term->div->n_row;
1783         case isl_dim_all:       return isl_dim_total(term->dim) + term->div->n_row;
1784         }
1785 }
1786
1787 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
1788 {
1789         return term ? term->dim->ctx : NULL;
1790 }
1791
1792 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
1793 {
1794         if (!term)
1795                 return;
1796         isl_int_set(*n, term->n);
1797 }
1798
1799 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
1800 {
1801         if (!term)
1802                 return;
1803         isl_int_set(*d, term->d);
1804 }
1805
1806 int isl_term_get_exp(__isl_keep isl_term *term,
1807         enum isl_dim_type type, unsigned pos)
1808 {
1809         if (!term)
1810                 return -1;
1811
1812         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
1813
1814         if (type >= isl_dim_set)
1815                 pos += isl_dim_size(term->dim, isl_dim_param);
1816         if (type >= isl_dim_div)
1817                 pos += isl_dim_size(term->dim, isl_dim_set);
1818
1819         return term->pow[pos];
1820 }
1821
1822 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
1823 {
1824         isl_basic_map *bmap;
1825         unsigned total;
1826         int k;
1827
1828         if (!term)
1829                 return NULL;
1830
1831         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
1832                         return NULL);
1833
1834         total = term->div->n_col - term->div->n_row - 2;
1835         /* No nested divs for now */
1836         isl_assert(term->dim->ctx,
1837                 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
1838                                         term->div->n_row) == -1,
1839                 return NULL);
1840
1841         bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
1842         if ((k = isl_basic_map_alloc_div(bmap)) < 0)
1843                 goto error;
1844
1845         isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
1846
1847         return isl_basic_map_div(bmap, k);
1848 error:
1849         isl_basic_map_free(bmap);
1850         return NULL;
1851 }
1852
1853 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
1854         int (*fn)(__isl_take isl_term *term, void *user),
1855         __isl_take isl_term *term, void *user)
1856 {
1857         int i;
1858         struct isl_upoly_rec *rec;
1859
1860         if (!up || !term)
1861                 goto error;
1862
1863         if (isl_upoly_is_zero(up))
1864                 return term;
1865
1866         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
1867         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
1868         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
1869
1870         if (isl_upoly_is_cst(up)) {
1871                 struct isl_upoly_cst *cst;
1872                 cst = isl_upoly_as_cst(up);
1873                 if (!cst)
1874                         goto error;
1875                 term = isl_term_cow(term);
1876                 if (!term)
1877                         goto error;
1878                 isl_int_set(term->n, cst->n);
1879                 isl_int_set(term->d, cst->d);
1880                 if (fn(isl_term_copy(term), user) < 0)
1881                         goto error;
1882                 return term;
1883         }
1884
1885         rec = isl_upoly_as_rec(up);
1886         if (!rec)
1887                 goto error;
1888
1889         for (i = 0; i < rec->n; ++i) {
1890                 term = isl_term_cow(term);
1891                 if (!term)
1892                         goto error;
1893                 term->pow[up->var] = i;
1894                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
1895                 if (!term)
1896                         goto error;
1897         }
1898         term->pow[up->var] = 0;
1899
1900         return term;
1901 error:
1902         isl_term_free(term);
1903         return NULL;
1904 }
1905
1906 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
1907         int (*fn)(__isl_take isl_term *term, void *user), void *user)
1908 {
1909         isl_term *term;
1910
1911         if (!qp)
1912                 return -1;
1913
1914         term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
1915         if (!term)
1916                 return -1;
1917
1918         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
1919
1920         isl_term_free(term);
1921
1922         return term ? 0 : -1;
1923 }
1924
1925 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
1926         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
1927                     void *user), void *user)
1928 {
1929         int i;
1930
1931         if (!pwqp)
1932                 return -1;
1933
1934         for (i = 0; i < pwqp->n; ++i)
1935                 if (fn(isl_set_copy(pwqp->p[i].set),
1936                                 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
1937                         return -1;
1938
1939         return 0;
1940 }
1941
1942 static int any_divs(__isl_keep isl_set *set)
1943 {
1944         int i;
1945
1946         if (!set)
1947                 return -1;
1948
1949         for (i = 0; i < set->n; ++i)
1950                 if (set->p[i]->n_div > 0)
1951                         return 1;
1952
1953         return 0;
1954 }
1955
1956 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
1957         __isl_take isl_dim *dim)
1958 {
1959         if (!qp || !dim)
1960                 goto error;
1961
1962         if (isl_dim_equal(qp->dim, dim)) {
1963                 isl_dim_free(dim);
1964                 return qp;
1965         }
1966
1967         qp = isl_qpolynomial_cow(qp);
1968         if (!qp)
1969                 goto error;
1970
1971         if (qp->div->n_row) {
1972                 int i;
1973                 int extra;
1974                 unsigned total;
1975                 int *exp;
1976
1977                 extra = isl_dim_size(dim, isl_dim_set) -
1978                                 isl_dim_size(qp->dim, isl_dim_set);
1979                 total = isl_dim_total(qp->dim);
1980                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1981                 if (!exp)
1982                         goto error;
1983                 for (i = 0; i < qp->div->n_row; ++i)
1984                         exp[i] = extra + i;
1985                 qp->upoly = expand(qp->upoly, exp, total);
1986                 free(exp);
1987                 if (!qp->upoly)
1988                         goto error;
1989                 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
1990                 if (!qp->div)
1991                         goto error;
1992                 for (i = 0; i < qp->div->n_row; ++i)
1993                         isl_seq_clr(qp->div->row[i] + 2 + total, extra);
1994         }
1995
1996         isl_dim_free(qp->dim);
1997         qp->dim = dim;
1998
1999         return qp;
2000 error:
2001         isl_dim_free(dim);
2002         isl_qpolynomial_free(qp);
2003         return NULL;
2004 }
2005
2006 static int foreach_lifted_subset(__isl_take isl_set *set,
2007         __isl_take isl_qpolynomial *qp,
2008         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2009                     void *user), void *user)
2010 {
2011         int i;
2012
2013         if (!set || !qp)
2014                 goto error;
2015
2016         for (i = 0; i < set->n; ++i) {
2017                 isl_set *lift;
2018                 isl_qpolynomial *copy;
2019
2020                 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2021                 lift = isl_set_lift(lift);
2022
2023                 copy = isl_qpolynomial_copy(qp);
2024                 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2025
2026                 if (fn(lift, copy, user) < 0)
2027                         goto error;
2028         }
2029
2030         isl_set_free(set);
2031         isl_qpolynomial_free(qp);
2032
2033         return 0;
2034 error:
2035         isl_set_free(set);
2036         isl_qpolynomial_free(qp);
2037         return -1;
2038 }
2039
2040 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2041         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2042                     void *user), void *user)
2043 {
2044         int i;
2045
2046         if (!pwqp)
2047                 return -1;
2048
2049         for (i = 0; i < pwqp->n; ++i) {
2050                 isl_set *set;
2051                 isl_qpolynomial *qp;
2052
2053                 set = isl_set_copy(pwqp->p[i].set);
2054                 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2055                 if (!any_divs(set)) {
2056                         if (fn(set, qp, user) < 0)
2057                                 return -1;
2058                         continue;
2059                 }
2060                 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2061                         return -1;
2062         }
2063
2064         return 0;
2065 }
2066
2067 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
2068         enum isl_fold type, __isl_take isl_dim *dim, int n)
2069 {
2070         isl_qpolynomial_fold *fold;
2071
2072         if (!dim)
2073                 goto error;
2074
2075         isl_assert(dim->ctx, n >= 0, goto error);
2076         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
2077                         sizeof(struct isl_qpolynomial_fold) +
2078                         (n - 1) * sizeof(struct isl_qpolynomial *));
2079         if (!fold)
2080                 goto error;
2081
2082         fold->ref = 1;
2083         fold->size = n;
2084         fold->n = 0;
2085         fold->type = type;
2086         fold->dim = dim;
2087
2088         return fold;
2089 error:
2090         isl_dim_free(dim);
2091         return NULL;
2092 }
2093
2094 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
2095         __isl_take isl_dim *dim)
2096 {
2097         return qpolynomial_fold_alloc(type, dim, 0);
2098 }
2099
2100 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
2101         enum isl_fold type, __isl_take isl_qpolynomial *qp)
2102 {
2103         isl_qpolynomial_fold *fold;
2104
2105         if (!qp)
2106                 return NULL;
2107
2108         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
2109         if (!fold)
2110                 goto error;
2111
2112         fold->qp[0] = qp;
2113         fold->n++;
2114
2115         return fold;
2116 error:
2117         isl_qpolynomial_fold_free(fold);
2118         isl_qpolynomial_free(qp);
2119         return NULL;
2120 }
2121
2122 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
2123         __isl_keep isl_qpolynomial_fold *fold)
2124 {
2125         if (!fold)
2126                 return NULL;
2127
2128         fold->ref++;
2129         return fold;
2130 }
2131
2132 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
2133 {
2134         int i;
2135
2136         if (!fold)
2137                 return;
2138         if (--fold->ref > 0)
2139                 return;
2140
2141         for (i = 0; i < fold->n; ++i)
2142                 isl_qpolynomial_free(fold->qp[i]);
2143         isl_dim_free(fold->dim);
2144         free(fold);
2145 }
2146
2147 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
2148 {
2149         if (!fold)
2150                 return -1;
2151
2152         return fold->n == 0;
2153 }
2154
2155 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
2156         __isl_take isl_qpolynomial_fold *fold1,
2157         __isl_take isl_qpolynomial_fold *fold2)
2158 {
2159         int i;
2160         struct isl_qpolynomial_fold *res = NULL;
2161
2162         if (!fold1 || !fold2)
2163                 goto error;
2164
2165         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
2166         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
2167                         goto error);
2168
2169         if (isl_qpolynomial_fold_is_empty(fold1)) {
2170                 isl_qpolynomial_fold_free(fold1);
2171                 return fold2;
2172         }
2173
2174         if (isl_qpolynomial_fold_is_empty(fold2)) {
2175                 isl_qpolynomial_fold_free(fold2);
2176                 return fold1;
2177         }
2178
2179         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
2180                                         fold1->n + fold2->n);
2181         if (!res)
2182                 goto error;
2183
2184         for (i = 0; i < fold1->n; ++i) {
2185                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
2186                 if (!res->qp[res->n])
2187                         goto error;
2188                 res->n++;
2189         }
2190
2191         for (i = 0; i < fold2->n; ++i) {
2192                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
2193                 if (!res->qp[res->n])
2194                         goto error;
2195                 res->n++;
2196         }
2197
2198         isl_qpolynomial_fold_free(fold1);
2199         isl_qpolynomial_fold_free(fold2);
2200
2201         return res;
2202 error:
2203         isl_qpolynomial_fold_free(res);
2204         isl_qpolynomial_fold_free(fold1);
2205         isl_qpolynomial_fold_free(fold2);
2206         return NULL;
2207 }
2208
2209 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2210         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
2211 {
2212         int i;
2213         isl_pw_qpolynomial_fold *pwf;
2214
2215         if (!pwqp)
2216                 return NULL;
2217
2218         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
2219
2220         for (i = 0; i < pwqp->n; ++i)
2221                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
2222                         isl_set_copy(pwqp->p[i].set),
2223                         isl_qpolynomial_fold_alloc(type,
2224                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
2225
2226         isl_pw_qpolynomial_free(pwqp);
2227
2228         return pwf;
2229 }