isl_polynomial.c

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