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 static __isl_give isl_pw_qpolynomial *pw_qpolynomial_alloc(__isl_take isl_dim *dim,
1283         int n)
1284 {
1285         struct isl_pw_qpolynomial *pwqp;
1286
1287         if (!dim)
1288                 return NULL;
1289         isl_assert(dim->ctx, n >= 0, return NULL);
1290         pwqp = isl_alloc(dim->ctx, struct isl_pw_qpolynomial,
1291                         sizeof(struct isl_pw_qpolynomial) +
1292                         (n - 1) * sizeof(struct isl_pw_qpolynomial_piece));
1293         if (!pwqp)
1294                 goto error;
1295
1296         pwqp->ref = 1;
1297         pwqp->size = n;
1298         pwqp->n = 0;
1299         pwqp->dim = dim;
1300         return pwqp;
1301 error:
1302         isl_dim_free(dim);
1303         return NULL;
1304 }
1305
1306 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_piece(
1307         __isl_take isl_pw_qpolynomial *pwqp,
1308         __isl_take isl_set *set, __isl_take isl_qpolynomial *qp)
1309 {
1310         if (!pwqp || !set || !qp)
1311                 goto error;
1312
1313         if (isl_set_fast_is_empty(set) || isl_qpolynomial_is_zero(qp)) {
1314                 isl_set_free(set);
1315                 isl_qpolynomial_free(qp);
1316                 return pwqp;
1317         }
1318
1319         isl_assert(set->ctx, isl_dim_equal(pwqp->dim, qp->dim), goto error);
1320         isl_assert(set->ctx, pwqp->n < pwqp->size, goto error);
1321
1322         pwqp->p[pwqp->n].set = set;
1323         pwqp->p[pwqp->n].qp = qp;
1324         pwqp->n++;
1325
1326         return pwqp;
1327 error:
1328         isl_pw_qpolynomial_free(pwqp);
1329         isl_set_free(set);
1330         isl_qpolynomial_free(qp);
1331         return NULL;
1332 }
1333
1334 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_alloc(__isl_take isl_set *set,
1335         __isl_take isl_qpolynomial *qp)
1336 {
1337         struct isl_pw_qpolynomial *pwqp;
1338
1339         if (!set || !qp)
1340                 goto error;
1341
1342         pwqp = pw_qpolynomial_alloc(isl_set_get_dim(set), 1);
1343
1344         return isl_pw_qpolynomial_add_piece(pwqp, set, qp);
1345 error:
1346         isl_set_free(set);
1347         isl_qpolynomial_free(qp);
1348         return NULL;
1349 }
1350
1351 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_zero(__isl_take isl_dim *dim)
1352 {
1353         return pw_qpolynomial_alloc(dim, 0);
1354 }
1355
1356 int isl_pw_qpolynomial_is_zero(__isl_keep isl_pw_qpolynomial *pwqp)
1357 {
1358         if (!pwqp)
1359                 return -1;
1360
1361         return pwqp->n == 0;
1362 }
1363
1364 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1365 {
1366         if (!pwqp)
1367                 return -1;
1368
1369         if (pwqp->n != -1)
1370                 return 0;
1371
1372         if (!isl_set_fast_is_universe(pwqp->p[0].set))
1373                 return 0;
1374
1375         return isl_qpolynomial_is_one(pwqp->p[0].qp);
1376 }
1377
1378 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_copy(
1379         __isl_keep isl_pw_qpolynomial *pwqp)
1380 {
1381         if (!pwqp)
1382                 return;
1383
1384         pwqp->ref++;
1385         return pwqp;
1386 }
1387
1388 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_dup(
1389         __isl_keep isl_pw_qpolynomial *pwqp)
1390 {
1391         int i;
1392         struct isl_pw_qpolynomial *dup;
1393
1394         if (!pwqp)
1395                 return NULL;
1396
1397         dup = pw_qpolynomial_alloc(isl_dim_copy(pwqp->dim), pwqp->n);
1398         if (!dup)
1399                 return NULL;
1400
1401         for (i = 0; i < pwqp->n; ++i)
1402                 dup = isl_pw_qpolynomial_add_piece(dup,
1403                                             isl_set_copy(pwqp->p[i].set),
1404                                             isl_qpolynomial_copy(pwqp->p[i].qp));
1405
1406         return dup;
1407 error:
1408         isl_pw_qpolynomial_free(dup);
1409         return NULL;
1410 }
1411
1412 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_cow(
1413         __isl_take isl_pw_qpolynomial *pwqp)
1414 {
1415         if (!pwqp)
1416                 return NULL;
1417
1418         if (pwqp->ref == 1)
1419                 return pwqp;
1420         pwqp->ref--;
1421         return isl_pw_qpolynomial_dup(pwqp);
1422 }
1423
1424 void isl_pw_qpolynomial_free(__isl_take isl_pw_qpolynomial *pwqp)
1425 {
1426         int i;
1427
1428         if (!pwqp)
1429                 return;
1430         if (--pwqp->ref > 0)
1431                 return;
1432
1433         for (i = 0; i < pwqp->n; ++i) {
1434                 isl_set_free(pwqp->p[i].set);
1435                 isl_qpolynomial_free(pwqp->p[i].qp);
1436         }
1437         isl_dim_free(pwqp->dim);
1438         free(pwqp);
1439 }
1440
1441 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
1442         __isl_take isl_pw_qpolynomial *pwqp1,
1443         __isl_take isl_pw_qpolynomial *pwqp2)
1444 {
1445         int i, j, n;
1446         struct isl_pw_qpolynomial *res;
1447         isl_set *set;
1448
1449         if (!pwqp1 || !pwqp2)
1450                 goto error;
1451
1452         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1453                         goto error);
1454
1455         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1456                 isl_pw_qpolynomial_free(pwqp1);
1457                 return pwqp2;
1458         }
1459
1460         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1461                 isl_pw_qpolynomial_free(pwqp2);
1462                 return pwqp1;
1463         }
1464
1465         n = (pwqp1->n + 1) * (pwqp2->n + 1);
1466         res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), n);
1467
1468         for (i = 0; i < pwqp1->n; ++i) {
1469                 set = isl_set_copy(pwqp1->p[i].set);
1470                 for (j = 0; j < pwqp2->n; ++j) {
1471                         struct isl_set *common;
1472                         struct isl_qpolynomial *sum;
1473                         set = isl_set_subtract(set,
1474                                         isl_set_copy(pwqp2->p[j].set));
1475                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1476                                                 isl_set_copy(pwqp2->p[j].set));
1477                         if (isl_set_fast_is_empty(common)) {
1478                                 isl_set_free(common);
1479                                 continue;
1480                         }
1481
1482                         sum = isl_qpolynomial_add(
1483                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1484                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1485
1486                         res = isl_pw_qpolynomial_add_piece(res, common, sum);
1487                 }
1488                 res = isl_pw_qpolynomial_add_piece(res, set,
1489                         isl_qpolynomial_copy(pwqp1->p[i].qp));
1490         }
1491
1492         for (j = 0; j < pwqp2->n; ++j) {
1493                 set = isl_set_copy(pwqp2->p[j].set);
1494                 for (i = 0; i < pwqp1->n; ++i)
1495                         set = isl_set_subtract(set,
1496                                         isl_set_copy(pwqp1->p[i].set));
1497                 res = isl_pw_qpolynomial_add_piece(res, set,
1498                         isl_qpolynomial_copy(pwqp2->p[j].qp));
1499         }
1500
1501         isl_pw_qpolynomial_free(pwqp1);
1502         isl_pw_qpolynomial_free(pwqp2);
1503
1504         return res;
1505 error:
1506         isl_pw_qpolynomial_free(pwqp1);
1507         isl_pw_qpolynomial_free(pwqp2);
1508         return NULL;
1509 }
1510
1511 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_disjoint(
1512         __isl_take isl_pw_qpolynomial *pwqp1,
1513         __isl_take isl_pw_qpolynomial *pwqp2)
1514 {
1515         int i;
1516         struct isl_pw_qpolynomial *res;
1517
1518         if (!pwqp1 || !pwqp2)
1519                 goto error;
1520
1521         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1522                         goto error);
1523
1524         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1525                 isl_pw_qpolynomial_free(pwqp1);
1526                 return pwqp2;
1527         }
1528
1529         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1530                 isl_pw_qpolynomial_free(pwqp2);
1531                 return pwqp1;
1532         }
1533
1534         res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), pwqp1->n + pwqp2->n);
1535
1536         for (i = 0; i < pwqp1->n; ++i)
1537                 res = isl_pw_qpolynomial_add_piece(res,
1538                                 isl_set_copy(pwqp1->p[i].set),
1539                                 isl_qpolynomial_copy(pwqp1->p[i].qp));
1540
1541         for (i = 0; i < pwqp2->n; ++i)
1542                 res = isl_pw_qpolynomial_add_piece(res,
1543                                 isl_set_copy(pwqp2->p[i].set),
1544                                 isl_qpolynomial_copy(pwqp2->p[i].qp));
1545
1546         isl_pw_qpolynomial_free(pwqp1);
1547         isl_pw_qpolynomial_free(pwqp2);
1548
1549         return res;
1550 error:
1551         isl_pw_qpolynomial_free(pwqp1);
1552         isl_pw_qpolynomial_free(pwqp2);
1553         return NULL;
1554 }
1555
1556 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1557         __isl_take isl_pw_qpolynomial *pwqp1,
1558         __isl_take isl_pw_qpolynomial *pwqp2)
1559 {
1560         int i, j, n;
1561         struct isl_pw_qpolynomial *res;
1562         isl_set *set;
1563
1564         if (!pwqp1 || !pwqp2)
1565                 goto error;
1566
1567         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1568                         goto error);
1569
1570         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1571                 isl_pw_qpolynomial_free(pwqp2);
1572                 return pwqp1;
1573         }
1574
1575         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1576                 isl_pw_qpolynomial_free(pwqp1);
1577                 return pwqp2;
1578         }
1579
1580         if (isl_pw_qpolynomial_is_one(pwqp1)) {
1581                 isl_pw_qpolynomial_free(pwqp1);
1582                 return pwqp2;
1583         }
1584
1585         if (isl_pw_qpolynomial_is_one(pwqp2)) {
1586                 isl_pw_qpolynomial_free(pwqp2);
1587                 return pwqp1;
1588         }
1589
1590         n = pwqp1->n * pwqp2->n;
1591         res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), n);
1592
1593         for (i = 0; i < pwqp1->n; ++i) {
1594                 for (j = 0; j < pwqp2->n; ++j) {
1595                         struct isl_set *common;
1596                         struct isl_qpolynomial *prod;
1597                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1598                                                 isl_set_copy(pwqp2->p[j].set));
1599                         if (isl_set_fast_is_empty(common)) {
1600                                 isl_set_free(common);
1601                                 continue;
1602                         }
1603
1604                         prod = isl_qpolynomial_mul(
1605                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1606                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1607
1608                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
1609                 }
1610         }
1611
1612         isl_pw_qpolynomial_free(pwqp1);
1613         isl_pw_qpolynomial_free(pwqp2);
1614
1615         return res;
1616 error:
1617         isl_pw_qpolynomial_free(pwqp1);
1618         isl_pw_qpolynomial_free(pwqp2);
1619         return NULL;
1620 }
1621
1622 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1623         __isl_take isl_pw_qpolynomial *pwqp)
1624 {
1625         int i, j, n;
1626         struct isl_pw_qpolynomial *res;
1627         isl_set *set;
1628
1629         if (!pwqp)
1630                 return NULL;
1631
1632         if (isl_pw_qpolynomial_is_zero(pwqp))
1633                 return pwqp;
1634
1635         pwqp = isl_pw_qpolynomial_cow(pwqp);
1636
1637         for (i = 0; i < pwqp->n; ++i) {
1638                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1639                 if (!pwqp->p[i].qp)
1640                         goto error;
1641         }
1642
1643         return pwqp;
1644 error:
1645         isl_pw_qpolynomial_free(pwqp);
1646         return NULL;
1647 }
1648
1649 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1650         __isl_take isl_pw_qpolynomial *pwqp1,
1651         __isl_take isl_pw_qpolynomial *pwqp2)
1652 {
1653         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1654 }
1655
1656 __isl_give struct isl_upoly *isl_upoly_eval(
1657         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1658 {
1659         int i;
1660         struct isl_upoly_rec *rec;
1661         struct isl_upoly *res;
1662         struct isl_upoly *base;
1663
1664         if (isl_upoly_is_cst(up)) {
1665                 isl_vec_free(vec);
1666                 return up;
1667         }
1668
1669         rec = isl_upoly_as_rec(up);
1670         if (!rec)
1671                 goto error;
1672
1673         isl_assert(up->ctx, rec->n >= 1, goto error);
1674
1675         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1676
1677         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1678                                 isl_vec_copy(vec));
1679
1680         for (i = rec->n - 2; i >= 0; --i) {
1681                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1682                 res = isl_upoly_sum(res,
1683                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1684                                                             isl_vec_copy(vec)));
1685         }
1686
1687         isl_upoly_free(base);
1688         isl_upoly_free(up);
1689         isl_vec_free(vec);
1690         return res;
1691 error:
1692         isl_upoly_free(up);
1693         isl_vec_free(vec);
1694         return NULL;
1695 }
1696
1697 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1698         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1699 {
1700         isl_vec *ext;
1701         struct isl_upoly *up;
1702
1703         if (!qp || !pnt)
1704                 goto error;
1705         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1706
1707         if (qp->div->n_row == 0)
1708                 ext = isl_vec_copy(pnt->vec);
1709         else {
1710                 int i;
1711                 unsigned dim = isl_dim_total(qp->dim);
1712                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1713                 if (!ext)
1714                         goto error;
1715
1716                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1717                 for (i = 0; i < qp->div->n_row; ++i) {
1718                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1719                                                 1 + dim + i, &ext->el[1+dim+i]);
1720                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1721                                         qp->div->row[i][0]);
1722                 }
1723         }
1724
1725         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1726         if (!up)
1727                 goto error;
1728
1729         isl_qpolynomial_free(qp);
1730         isl_point_free(pnt);
1731
1732         qp = isl_qpolynomial_alloc(isl_dim_set_alloc(up->ctx, 0, 0), 0);
1733         if (!qp)
1734                 isl_upoly_free(up);
1735         else
1736                 qp->upoly = up;
1737
1738         return qp;
1739 error:
1740         isl_qpolynomial_free(qp);
1741         isl_point_free(pnt);
1742         return NULL;
1743 }
1744
1745 __isl_give isl_qpolynomial *isl_pw_qpolynomial_eval(
1746         __isl_take isl_pw_qpolynomial *pwqp, __isl_take isl_point *pnt)
1747 {
1748         int i;
1749         int found;
1750         isl_qpolynomial *qp;
1751
1752         if (!pwqp || !pnt)
1753                 goto error;
1754         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, pwqp->dim), goto error);
1755
1756         for (i = 0; i < pwqp->n; ++i) {
1757                 found = isl_set_contains_point(pwqp->p[i].set, pnt);
1758                 if (found < 0)
1759                         goto error;
1760                 if (found)
1761                         break;
1762         }
1763         if (found)
1764                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(pwqp->p[i].qp),
1765                                             isl_point_copy(pnt));
1766         else
1767                 qp = isl_qpolynomial_zero(isl_dim_copy(pwqp->dim));
1768         isl_pw_qpolynomial_free(pwqp);
1769         isl_point_free(pnt);
1770         return qp;
1771 error:
1772         isl_pw_qpolynomial_free(pwqp);
1773         isl_point_free(pnt);
1774         return NULL;
1775 }
1776
1777 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1778         enum isl_dim_type dst_type, unsigned dst_pos,
1779         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1780 {
1781         if (!qp)
1782                 return NULL;
1783
1784         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1785         if (!qp->dim)
1786                 goto error;
1787
1788         /* Do something to polynomials when needed; later */
1789
1790         return qp;
1791 error:
1792         isl_qpolynomial_free(qp);
1793         return NULL;
1794 }
1795
1796 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1797         __isl_take isl_pw_qpolynomial *pwqp,
1798         enum isl_dim_type dst_type, unsigned dst_pos,
1799         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1800 {
1801         int i;
1802
1803         if (!pwqp)
1804                 return NULL;
1805
1806         pwqp->dim = isl_dim_move(pwqp->dim,
1807                                     dst_type, dst_pos, src_type, src_pos, n);
1808         if (!pwqp->dim)
1809                 goto error;
1810
1811         for (i = 0; i < pwqp->n; ++i) {
1812                 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1813                                                 src_type, src_pos, n);
1814                 if (!pwqp->p[i].set)
1815                         goto error;
1816                 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1817                                         dst_type, dst_pos, src_type, src_pos, n);
1818                 if (!pwqp->p[i].qp)
1819                         goto error;
1820         }
1821
1822         return pwqp;
1823 error:
1824         isl_pw_qpolynomial_free(pwqp);
1825         return NULL;
1826 }