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