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