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