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, *at = 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         at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1089         reordering = isl_alloc_array(qp->div->ctx, int, qp->div->n_col - 2);
1090         if (!array || !pos || !at || !reordering)
1091                 goto error;
1092
1093         for (i = 0; i < qp->div->n_row; ++i) {
1094                 array[i].div = qp->div;
1095                 array[i].row = i;
1096                 pos[i] = i;
1097                 at[i] = i;
1098         }
1099
1100         qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1101                 div_sort_cmp);
1102
1103         for (i = 0; i < div_pos; ++i)
1104                 reordering[i] = i;
1105
1106         for (i = 0; i < qp->div->n_row; ++i)
1107                 reordering[div_pos + array[i].row] = div_pos + i;
1108
1109         for (i = 0; i < qp->div->n_row; ++i) {
1110                 if (pos[array[i].row] == i)
1111                         continue;
1112                 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1113                 pos[at[i]] = pos[array[i].row];
1114                 at[pos[array[i].row]] = at[i];
1115                 at[i] = array[i].row;
1116                 pos[array[i].row] = i;
1117         }
1118
1119         qp->upoly = reorder(qp->upoly, reordering);
1120
1121         if (!qp->upoly || !qp->div)
1122                 goto error;
1123
1124         free(at);
1125         free(pos);
1126         free(array);
1127         free(reordering);
1128
1129         return qp;
1130 error:
1131         free(at);
1132         free(pos);
1133         free(array);
1134         free(reordering);
1135         isl_qpolynomial_free(qp);
1136         return NULL;
1137 }
1138
1139 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
1140         __isl_keep isl_mat *div2, int *exp1, int *exp2)
1141 {
1142         int i, j, k;
1143         isl_mat *div = NULL;
1144         unsigned d = div1->n_col - div1->n_row;
1145
1146         div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
1147                                 d + div1->n_row + div2->n_row);
1148         if (!div)
1149                 return NULL;
1150
1151         for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
1152                 int cmp;
1153
1154                 expand_row(div, k, div1, i, exp1);
1155                 expand_row(div, k + 1, div2, j, exp2);
1156
1157                 cmp = cmp_row(div, k, k + 1);
1158                 if (cmp == 0) {
1159                         exp1[i++] = k;
1160                         exp2[j++] = k;
1161                 } else if (cmp < 0) {
1162                         exp1[i++] = k;
1163                 } else {
1164                         exp2[j++] = k;
1165                         isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
1166                 }
1167         }
1168         for (; i < div1->n_row; ++i, ++k) {
1169                 expand_row(div, k, div1, i, exp1);
1170                 exp1[i] = k;
1171         }
1172         for (; j < div2->n_row; ++j, ++k) {
1173                 expand_row(div, k, div2, j, exp2);
1174                 exp2[j] = k;
1175         }
1176
1177         div->n_row = k;
1178         div->n_col = d + k;
1179
1180         return div;
1181 }
1182
1183 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1184         int *exp, int first)
1185 {
1186         int i;
1187         struct isl_upoly_rec *rec;
1188
1189         if (isl_upoly_is_cst(up))
1190                 return up;
1191
1192         if (up->var < first)
1193                 return up;
1194
1195         if (exp[up->var - first] == up->var - first)
1196                 return up;
1197
1198         up = isl_upoly_cow(up);
1199         if (!up)
1200                 goto error;
1201
1202         up->var = exp[up->var - first] + first;
1203
1204         rec = isl_upoly_as_rec(up);
1205         if (!rec)
1206                 goto error;
1207
1208         for (i = 0; i < rec->n; ++i) {
1209                 rec->p[i] = expand(rec->p[i], exp, first);
1210                 if (!rec->p[i])
1211                         goto error;
1212         }
1213
1214         return up;
1215 error:
1216         isl_upoly_free(up);
1217         return NULL;
1218 }
1219
1220 static __isl_give isl_qpolynomial *with_merged_divs(
1221         __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1222                                           __isl_take isl_qpolynomial *qp2),
1223         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1224 {
1225         int *exp1 = NULL;
1226         int *exp2 = NULL;
1227         isl_mat *div = NULL;
1228
1229         qp1 = isl_qpolynomial_cow(qp1);
1230         qp2 = isl_qpolynomial_cow(qp2);
1231
1232         if (!qp1 || !qp2)
1233                 goto error;
1234
1235         isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1236                                 qp1->div->n_col >= qp2->div->n_col, goto error);
1237
1238         exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1239         exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1240         if (!exp1 || !exp2)
1241                 goto error;
1242
1243         div = merge_divs(qp1->div, qp2->div, exp1, exp2);
1244         if (!div)
1245                 goto error;
1246
1247         isl_mat_free(qp1->div);
1248         qp1->div = isl_mat_copy(div);
1249         isl_mat_free(qp2->div);
1250         qp2->div = isl_mat_copy(div);
1251
1252         qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1253         qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1254
1255         if (!qp1->upoly || !qp2->upoly)
1256                 goto error;
1257
1258         isl_mat_free(div);
1259         free(exp1);
1260         free(exp2);
1261
1262         return fn(qp1, qp2);
1263 error:
1264         isl_mat_free(div);
1265         free(exp1);
1266         free(exp2);
1267         isl_qpolynomial_free(qp1);
1268         isl_qpolynomial_free(qp2);
1269         return NULL;
1270 }
1271
1272 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1273         __isl_take isl_qpolynomial *qp2)
1274 {
1275         qp1 = isl_qpolynomial_cow(qp1);
1276
1277         if (!qp1 || !qp2)
1278                 goto error;
1279
1280         if (qp1->div->n_row < qp2->div->n_row)
1281                 return isl_qpolynomial_add(qp2, qp1);
1282
1283         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1284         if (!compatible_divs(qp1->div, qp2->div))
1285                 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1286
1287         qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1288         if (!qp1->upoly)
1289                 goto error;
1290
1291         isl_qpolynomial_free(qp2);
1292
1293         return qp1;
1294 error:
1295         isl_qpolynomial_free(qp1);
1296         isl_qpolynomial_free(qp2);
1297         return NULL;
1298 }
1299
1300 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1301         __isl_keep isl_set *dom,
1302         __isl_take isl_qpolynomial *qp1,
1303         __isl_take isl_qpolynomial *qp2)
1304 {
1305         return isl_qpolynomial_add(qp1, qp2);
1306 }
1307
1308 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1309         __isl_take isl_qpolynomial *qp2)
1310 {
1311         return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1312 }
1313
1314 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1315 {
1316         qp = isl_qpolynomial_cow(qp);
1317
1318         if (!qp)
1319                 return NULL;
1320
1321         qp->upoly = isl_upoly_neg(qp->upoly);
1322         if (!qp->upoly)
1323                 goto error;
1324
1325         return qp;
1326 error:
1327         isl_qpolynomial_free(qp);
1328         return NULL;
1329 }
1330
1331 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1332         __isl_take isl_qpolynomial *qp2)
1333 {
1334         qp1 = isl_qpolynomial_cow(qp1);
1335
1336         if (!qp1 || !qp2)
1337                 goto error;
1338
1339         if (qp1->div->n_row < qp2->div->n_row)
1340                 return isl_qpolynomial_mul(qp2, qp1);
1341
1342         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1343         if (!compatible_divs(qp1->div, qp2->div))
1344                 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1345
1346         qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1347         if (!qp1->upoly)
1348                 goto error;
1349
1350         isl_qpolynomial_free(qp2);
1351
1352         return qp1;
1353 error:
1354         isl_qpolynomial_free(qp1);
1355         isl_qpolynomial_free(qp2);
1356         return NULL;
1357 }
1358
1359 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1360 {
1361         return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1362 }
1363
1364 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1365 {
1366         return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1367 }
1368
1369 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty(__isl_take isl_dim *dim)
1370 {
1371         return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1372 }
1373
1374 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1375 {
1376         return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1377 }
1378
1379 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1380         isl_int v)
1381 {
1382         struct isl_qpolynomial *qp;
1383         struct isl_upoly_cst *cst;
1384
1385         qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1386         if (!qp)
1387                 return NULL;
1388
1389         cst = isl_upoly_as_cst(qp->upoly);
1390         isl_int_set(cst->n, v);
1391
1392         return qp;
1393 }
1394
1395 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1396         isl_int *n, isl_int *d)
1397 {
1398         struct isl_upoly_cst *cst;
1399
1400         if (!qp)
1401                 return -1;
1402
1403         if (!isl_upoly_is_cst(qp->upoly))
1404                 return 0;
1405
1406         cst = isl_upoly_as_cst(qp->upoly);
1407         if (!cst)
1408                 return -1;
1409
1410         if (n)
1411                 isl_int_set(*n, cst->n);
1412         if (d)
1413                 isl_int_set(*d, cst->d);
1414
1415         return 1;
1416 }
1417
1418 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1419 {
1420         int is_cst;
1421         struct isl_upoly_rec *rec;
1422
1423         if (!up)
1424                 return -1;
1425
1426         if (up->var < 0)
1427                 return 1;
1428
1429         rec = isl_upoly_as_rec(up);
1430         if (!rec)
1431                 return -1;
1432
1433         if (rec->n > 2)
1434                 return 0;
1435
1436         isl_assert(up->ctx, rec->n > 1, return -1);
1437
1438         is_cst = isl_upoly_is_cst(rec->p[1]);
1439         if (is_cst < 0)
1440                 return -1;
1441         if (!is_cst)
1442                 return 0;
1443
1444         return isl_upoly_is_affine(rec->p[0]);
1445 }
1446
1447 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1448 {
1449         if (!qp)
1450                 return -1;
1451
1452         if (qp->div->n_row > 0)
1453                 return 0;
1454
1455         return isl_upoly_is_affine(qp->upoly);
1456 }
1457
1458 static void update_coeff(__isl_keep isl_vec *aff,
1459         __isl_keep struct isl_upoly_cst *cst, int pos)
1460 {
1461         isl_int gcd;
1462         isl_int f;
1463
1464         if (isl_int_is_zero(cst->n))
1465                 return;
1466
1467         isl_int_init(gcd);
1468         isl_int_init(f);
1469         isl_int_gcd(gcd, cst->d, aff->el[0]);
1470         isl_int_divexact(f, cst->d, gcd);
1471         isl_int_divexact(gcd, aff->el[0], gcd);
1472         isl_seq_scale(aff->el, aff->el, f, aff->size);
1473         isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1474         isl_int_clear(gcd);
1475         isl_int_clear(f);
1476 }
1477
1478 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1479         __isl_keep isl_vec *aff)
1480 {
1481         struct isl_upoly_cst *cst;
1482         struct isl_upoly_rec *rec;
1483
1484         if (!up || !aff)
1485                 return -1;
1486
1487         if (up->var < 0) {
1488                 struct isl_upoly_cst *cst;
1489
1490                 cst = isl_upoly_as_cst(up);
1491                 if (!cst)
1492                         return -1;
1493                 update_coeff(aff, cst, 0);
1494                 return 0;
1495         }
1496
1497         rec = isl_upoly_as_rec(up);
1498         if (!rec)
1499                 return -1;
1500         isl_assert(up->ctx, rec->n == 2, return -1);
1501
1502         cst = isl_upoly_as_cst(rec->p[1]);
1503         if (!cst)
1504                 return -1;
1505         update_coeff(aff, cst, 1 + up->var);
1506
1507         return isl_upoly_update_affine(rec->p[0], aff);
1508 }
1509
1510 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1511         __isl_keep isl_qpolynomial *qp)
1512 {
1513         isl_vec *aff;
1514         unsigned d;
1515
1516         if (!qp)
1517                 return NULL;
1518
1519         isl_assert(qp->div->ctx, qp->div->n_row == 0, return NULL);
1520         d = isl_dim_total(qp->dim);
1521         aff = isl_vec_alloc(qp->div->ctx, 2 + d);
1522         if (!aff)
1523                 return NULL;
1524
1525         isl_seq_clr(aff->el + 1, 1 + d);
1526         isl_int_set_si(aff->el[0], 1);
1527
1528         if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1529                 goto error;
1530
1531         return aff;
1532 error:
1533         isl_vec_free(aff);
1534         return NULL;
1535 }
1536
1537 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
1538         __isl_keep isl_qpolynomial *qp2)
1539 {
1540         if (!qp1 || !qp2)
1541                 return -1;
1542
1543         return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1544 }
1545
1546 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1547 {
1548         int i;
1549         struct isl_upoly_rec *rec;
1550
1551         if (isl_upoly_is_cst(up)) {
1552                 struct isl_upoly_cst *cst;
1553                 cst = isl_upoly_as_cst(up);
1554                 if (!cst)
1555                         return;
1556                 isl_int_lcm(*d, *d, cst->d);
1557                 return;
1558         }
1559
1560         rec = isl_upoly_as_rec(up);
1561         if (!rec)
1562                 return;
1563
1564         for (i = 0; i < rec->n; ++i)
1565                 upoly_update_den(rec->p[i], d);
1566 }
1567
1568 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1569 {
1570         isl_int_set_si(*d, 1);
1571         if (!qp)
1572                 return;
1573         upoly_update_den(qp->upoly, d);
1574 }
1575
1576 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1577         int pos, int power)
1578 {
1579         struct isl_ctx *ctx;
1580
1581         if (!dim)
1582                 return NULL;
1583
1584         ctx = dim->ctx;
1585
1586         return isl_qpolynomial_alloc(dim, 0, isl_upoly_pow(ctx, pos, power));
1587 }
1588
1589 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1590         enum isl_dim_type type, unsigned pos)
1591 {
1592         if (!dim)
1593                 return NULL;
1594
1595         isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1596         isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1597
1598         if (type == isl_dim_set)
1599                 pos += isl_dim_size(dim, isl_dim_param);
1600
1601         return isl_qpolynomial_pow(dim, pos, 1);
1602 error:
1603         isl_dim_free(dim);
1604         return NULL;
1605 }
1606
1607 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1608         int power)
1609 {
1610         struct isl_qpolynomial *qp = NULL;
1611         struct isl_upoly_rec *rec;
1612         struct isl_upoly_cst *cst;
1613         int i;
1614         int pos;
1615
1616         if (!div)
1617                 return NULL;
1618         isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1619
1620         pos = isl_dim_total(div->bmap->dim);
1621         rec = isl_upoly_alloc_rec(div->ctx, pos, 1 + power);
1622         qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1,
1623                                    &rec->up);
1624         if (!qp)
1625                 goto error;
1626
1627         isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1628         isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1629
1630         for (i = 0; i < 1 + power; ++i) {
1631                 rec->p[i] = isl_upoly_zero(div->ctx);
1632                 if (!rec->p[i])
1633                         goto error;
1634                 rec->n++;
1635         }
1636         cst = isl_upoly_as_cst(rec->p[power]);
1637         isl_int_set_si(cst->n, 1);
1638
1639         isl_div_free(div);
1640
1641         return qp;
1642 error:
1643         isl_qpolynomial_free(qp);
1644         isl_div_free(div);
1645         return NULL;
1646 }
1647
1648 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1649 {
1650         return isl_qpolynomial_div_pow(div, 1);
1651 }
1652
1653 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1654         const isl_int n, const isl_int d)
1655 {
1656         struct isl_qpolynomial *qp;
1657         struct isl_upoly_cst *cst;
1658
1659         qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1660         if (!qp)
1661                 return NULL;
1662
1663         cst = isl_upoly_as_cst(qp->upoly);
1664         isl_int_set(cst->n, n);
1665         isl_int_set(cst->d, d);
1666
1667         return qp;
1668 }
1669
1670 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
1671 {
1672         struct isl_upoly_rec *rec;
1673         int i;
1674
1675         if (!up)
1676                 return -1;
1677
1678         if (isl_upoly_is_cst(up))
1679                 return 0;
1680
1681         if (up->var < d)
1682                 active[up->var] = 1;
1683
1684         rec = isl_upoly_as_rec(up);
1685         for (i = 0; i < rec->n; ++i)
1686                 if (up_set_active(rec->p[i], active, d) < 0)
1687                         return -1;
1688
1689         return 0;
1690 }
1691
1692 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
1693 {
1694         int i, j;
1695         int d = isl_dim_total(qp->dim);
1696
1697         if (!qp || !active)
1698                 return -1;
1699
1700         for (i = 0; i < d; ++i)
1701                 for (j = 0; j < qp->div->n_row; ++j) {
1702                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
1703                                 continue;
1704                         active[i] = 1;
1705                         break;
1706                 }
1707
1708         return up_set_active(qp->upoly, active, d);
1709 }
1710
1711 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
1712         enum isl_dim_type type, unsigned first, unsigned n)
1713 {
1714         int i;
1715         int *active = NULL;
1716         int involves = 0;
1717
1718         if (!qp)
1719                 return -1;
1720         if (n == 0)
1721                 return 0;
1722
1723         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1724                         return -1);
1725         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1726                                  type == isl_dim_set, return -1);
1727
1728         active = isl_calloc_array(set->ctx, int, isl_dim_total(qp->dim));
1729         if (set_active(qp, active) < 0)
1730                 goto error;
1731
1732         if (type == isl_dim_set)
1733                 first += isl_dim_size(qp->dim, isl_dim_param);
1734         for (i = 0; i < n; ++i)
1735                 if (active[first + i]) {
1736                         involves = 1;
1737                         break;
1738                 }
1739
1740         free(active);
1741
1742         return involves;
1743 error:
1744         free(active);
1745         return -1;
1746 }
1747
1748 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
1749         unsigned first, unsigned n)
1750 {
1751         int i;
1752         struct isl_upoly_rec *rec;
1753
1754         if (!up)
1755                 return NULL;
1756         if (n == 0 || up->var < 0 || up->var < first)
1757                 return up;
1758         if (up->var < first + n) {
1759                 up = replace_by_constant_term(up);
1760                 return isl_upoly_drop(up, first, n);
1761         }
1762         up = isl_upoly_cow(up);
1763         if (!up)
1764                 return NULL;
1765         up->var -= n;
1766         rec = isl_upoly_as_rec(up);
1767         if (!rec)
1768                 goto error;
1769
1770         for (i = 0; i < rec->n; ++i) {
1771                 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
1772                 if (!rec->p[i])
1773                         goto error;
1774         }
1775
1776         return up;
1777 error:
1778         isl_upoly_free(up);
1779         return NULL;
1780 }
1781
1782 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
1783         __isl_take isl_qpolynomial *qp,
1784         enum isl_dim_type type, unsigned first, unsigned n)
1785 {
1786         if (!qp)
1787                 return NULL;
1788         if (n == 0 && !isl_dim_get_tuple_name(qp->dim, type))
1789                 return qp;
1790
1791         qp = isl_qpolynomial_cow(qp);
1792         if (!qp)
1793                 return NULL;
1794
1795         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1796                         goto error);
1797         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1798                                  type == isl_dim_set, goto error);
1799
1800         qp->dim = isl_dim_drop(qp->dim, type, first, n);
1801         if (!qp->dim)
1802                 goto error;
1803
1804         if (type == isl_dim_set)
1805                 first += isl_dim_size(qp->dim, isl_dim_param);
1806
1807         qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
1808         if (!qp->div)
1809                 goto error;
1810
1811         qp->upoly = isl_upoly_drop(qp->upoly, first, n);
1812         if (!qp->upoly)
1813                 goto error;
1814
1815         return qp;
1816 error:
1817         isl_qpolynomial_free(qp);
1818         return NULL;
1819 }
1820
1821 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1822         unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1823 {
1824         int i;
1825         struct isl_upoly_rec *rec;
1826         struct isl_upoly *base, *res;
1827
1828         if (!up)
1829                 return NULL;
1830
1831         if (isl_upoly_is_cst(up))
1832                 return up;
1833
1834         if (up->var < first)
1835                 return up;
1836
1837         rec = isl_upoly_as_rec(up);
1838         if (!rec)
1839                 goto error;
1840
1841         isl_assert(up->ctx, rec->n >= 1, goto error);
1842
1843         if (up->var >= first + n)
1844                 base = isl_upoly_pow(up->ctx, up->var, 1);
1845         else
1846                 base = isl_upoly_copy(subs[up->var - first]);
1847
1848         res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
1849         for (i = rec->n - 2; i >= 0; --i) {
1850                 struct isl_upoly *t;
1851                 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
1852                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1853                 res = isl_upoly_sum(res, t);
1854         }
1855
1856         isl_upoly_free(base);
1857         isl_upoly_free(up);
1858
1859         return res;
1860 error:
1861         isl_upoly_free(up);
1862         return NULL;
1863 }
1864
1865 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
1866         isl_int denom, unsigned len)
1867 {
1868         int i;
1869         struct isl_upoly *up;
1870
1871         isl_assert(ctx, len >= 1, return NULL);
1872
1873         up = isl_upoly_rat_cst(ctx, f[0], denom);
1874         for (i = 0; i < len - 1; ++i) {
1875                 struct isl_upoly *t;
1876                 struct isl_upoly *c;
1877
1878                 if (isl_int_is_zero(f[1 + i]))
1879                         continue;
1880
1881                 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
1882                 t = isl_upoly_pow(ctx, i, 1);
1883                 t = isl_upoly_mul(c, t);
1884                 up = isl_upoly_sum(up, t);
1885         }
1886
1887         return up;
1888 }
1889
1890 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
1891         __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
1892 {
1893         int i, j, k;
1894         isl_int denom;
1895         unsigned total;
1896         struct isl_upoly *up;
1897
1898         if (!eq)
1899                 goto error;
1900         if (eq->n_eq == 0) {
1901                 isl_basic_set_free(eq);
1902                 return qp;
1903         }
1904
1905         qp = isl_qpolynomial_cow(qp);
1906         if (!qp)
1907                 goto error;
1908         qp->div = isl_mat_cow(qp->div);
1909         if (!qp->div)
1910                 goto error;
1911
1912         total = 1 + isl_dim_total(eq->dim);
1913         isl_int_init(denom);
1914         for (i = 0; i < eq->n_eq; ++i) {
1915                 j = isl_seq_last_non_zero(eq->eq[i], total);
1916                 if (j < 0 || j == 0)
1917                         continue;
1918
1919                 for (k = 0; k < qp->div->n_row; ++k) {
1920                         if (isl_int_is_zero(qp->div->row[k][1 + j]))
1921                                 continue;
1922                         isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
1923                                         &qp->div->row[k][0]);
1924                         isl_seq_normalize(qp->div->ctx,
1925                                           qp->div->row[k], 1 + total);
1926                 }
1927
1928                 if (isl_int_is_pos(eq->eq[i][j]))
1929                         isl_seq_neg(eq->eq[i], eq->eq[i], total);
1930                 isl_int_abs(denom, eq->eq[i][j]);
1931                 isl_int_set_si(eq->eq[i][j], 0);
1932
1933                 up = isl_upoly_from_affine(qp->dim->ctx,
1934                                                    eq->eq[i], denom, total);
1935                 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
1936                 isl_upoly_free(up);
1937         }
1938         isl_int_clear(denom);
1939
1940         if (!qp->upoly)
1941                 goto error;
1942
1943         isl_basic_set_free(eq);
1944
1945         qp = sort_divs(qp);
1946
1947         return qp;
1948 error:
1949         isl_basic_set_free(eq);
1950         isl_qpolynomial_free(qp);
1951         return NULL;
1952 }
1953
1954 #undef PW
1955 #define PW isl_pw_qpolynomial
1956 #undef EL
1957 #define EL isl_qpolynomial
1958 #undef IS_ZERO
1959 #define IS_ZERO is_zero
1960 #undef FIELD
1961 #define FIELD qp
1962
1963 #include <isl_pw_templ.c>
1964
1965 #undef UNION
1966 #define UNION isl_union_pw_qpolynomial
1967 #undef PART
1968 #define PART isl_pw_qpolynomial
1969 #undef PARTS
1970 #define PARTS pw_qpolynomial
1971
1972 #include <isl_union_templ.c>
1973
1974 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1975 {
1976         if (!pwqp)
1977                 return -1;
1978
1979         if (pwqp->n != -1)
1980                 return 0;
1981
1982         if (!isl_set_fast_is_universe(pwqp->p[0].set))
1983                 return 0;
1984
1985         return isl_qpolynomial_is_one(pwqp->p[0].qp);
1986 }
1987
1988 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1989         __isl_take isl_pw_qpolynomial *pwqp1,
1990         __isl_take isl_pw_qpolynomial *pwqp2)
1991 {
1992         int i, j, n;
1993         struct isl_pw_qpolynomial *res;
1994         isl_set *set;
1995
1996         if (!pwqp1 || !pwqp2)
1997                 goto error;
1998
1999         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
2000                         goto error);
2001
2002         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2003                 isl_pw_qpolynomial_free(pwqp2);
2004                 return pwqp1;
2005         }
2006
2007         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2008                 isl_pw_qpolynomial_free(pwqp1);
2009                 return pwqp2;
2010         }
2011
2012         if (isl_pw_qpolynomial_is_one(pwqp1)) {
2013                 isl_pw_qpolynomial_free(pwqp1);
2014                 return pwqp2;
2015         }
2016
2017         if (isl_pw_qpolynomial_is_one(pwqp2)) {
2018                 isl_pw_qpolynomial_free(pwqp2);
2019                 return pwqp1;
2020         }
2021
2022         n = pwqp1->n * pwqp2->n;
2023         res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
2024
2025         for (i = 0; i < pwqp1->n; ++i) {
2026                 for (j = 0; j < pwqp2->n; ++j) {
2027                         struct isl_set *common;
2028                         struct isl_qpolynomial *prod;
2029                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2030                                                 isl_set_copy(pwqp2->p[j].set));
2031                         if (isl_set_fast_is_empty(common)) {
2032                                 isl_set_free(common);
2033                                 continue;
2034                         }
2035
2036                         prod = isl_qpolynomial_mul(
2037                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
2038                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
2039
2040                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
2041                 }
2042         }
2043
2044         isl_pw_qpolynomial_free(pwqp1);
2045         isl_pw_qpolynomial_free(pwqp2);
2046
2047         return res;
2048 error:
2049         isl_pw_qpolynomial_free(pwqp1);
2050         isl_pw_qpolynomial_free(pwqp2);
2051         return NULL;
2052 }
2053
2054 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
2055         __isl_take isl_pw_qpolynomial *pwqp)
2056 {
2057         int i;
2058
2059         if (!pwqp)
2060                 return NULL;
2061
2062         if (isl_pw_qpolynomial_is_zero(pwqp))
2063                 return pwqp;
2064
2065         pwqp = isl_pw_qpolynomial_cow(pwqp);
2066         if (!pwqp)
2067                 return NULL;
2068
2069         for (i = 0; i < pwqp->n; ++i) {
2070                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
2071                 if (!pwqp->p[i].qp)
2072                         goto error;
2073         }
2074
2075         return pwqp;
2076 error:
2077         isl_pw_qpolynomial_free(pwqp);
2078         return NULL;
2079 }
2080
2081 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
2082         __isl_take isl_pw_qpolynomial *pwqp1,
2083         __isl_take isl_pw_qpolynomial *pwqp2)
2084 {
2085         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
2086 }
2087
2088 __isl_give struct isl_upoly *isl_upoly_eval(
2089         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
2090 {
2091         int i;
2092         struct isl_upoly_rec *rec;
2093         struct isl_upoly *res;
2094         struct isl_upoly *base;
2095
2096         if (isl_upoly_is_cst(up)) {
2097                 isl_vec_free(vec);
2098                 return up;
2099         }
2100
2101         rec = isl_upoly_as_rec(up);
2102         if (!rec)
2103                 goto error;
2104
2105         isl_assert(up->ctx, rec->n >= 1, goto error);
2106
2107         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
2108
2109         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2110                                 isl_vec_copy(vec));
2111
2112         for (i = rec->n - 2; i >= 0; --i) {
2113                 res = isl_upoly_mul(res, isl_upoly_copy(base));
2114                 res = isl_upoly_sum(res,
2115                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2116                                                             isl_vec_copy(vec)));
2117         }
2118
2119         isl_upoly_free(base);
2120         isl_upoly_free(up);
2121         isl_vec_free(vec);
2122         return res;
2123 error:
2124         isl_upoly_free(up);
2125         isl_vec_free(vec);
2126         return NULL;
2127 }
2128
2129 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
2130         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
2131 {
2132         isl_vec *ext;
2133         struct isl_upoly *up;
2134         isl_dim *dim;
2135
2136         if (!qp || !pnt)
2137                 goto error;
2138         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
2139
2140         if (qp->div->n_row == 0)
2141                 ext = isl_vec_copy(pnt->vec);
2142         else {
2143                 int i;
2144                 unsigned dim = isl_dim_total(qp->dim);
2145                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2146                 if (!ext)
2147                         goto error;
2148
2149                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2150                 for (i = 0; i < qp->div->n_row; ++i) {
2151                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2152                                                 1 + dim + i, &ext->el[1+dim+i]);
2153                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2154                                         qp->div->row[i][0]);
2155                 }
2156         }
2157
2158         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2159         if (!up)
2160                 goto error;
2161
2162         dim = isl_dim_copy(qp->dim);
2163         isl_qpolynomial_free(qp);
2164         isl_point_free(pnt);
2165
2166         return isl_qpolynomial_alloc(dim, 0, up);
2167 error:
2168         isl_qpolynomial_free(qp);
2169         isl_point_free(pnt);
2170         return NULL;
2171 }
2172
2173 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2174         __isl_keep struct isl_upoly_cst *cst2)
2175 {
2176         int cmp;
2177         isl_int t;
2178         isl_int_init(t);
2179         isl_int_mul(t, cst1->n, cst2->d);
2180         isl_int_submul(t, cst2->n, cst1->d);
2181         cmp = isl_int_sgn(t);
2182         isl_int_clear(t);
2183         return cmp;
2184 }
2185
2186 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial *qp1,
2187         __isl_keep isl_qpolynomial *qp2)
2188 {
2189         struct isl_upoly_cst *cst1, *cst2;
2190
2191         if (!qp1 || !qp2)
2192                 return -1;
2193         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), return -1);
2194         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), return -1);
2195         if (isl_qpolynomial_is_nan(qp1))
2196                 return -1;
2197         if (isl_qpolynomial_is_nan(qp2))
2198                 return -1;
2199         cst1 = isl_upoly_as_cst(qp1->upoly);
2200         cst2 = isl_upoly_as_cst(qp2->upoly);
2201
2202         return isl_upoly_cmp(cst1, cst2) <= 0;
2203 }
2204
2205 __isl_give isl_qpolynomial *isl_qpolynomial_min_cst(
2206         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2207 {
2208         struct isl_upoly_cst *cst1, *cst2;
2209         int cmp;
2210
2211         if (!qp1 || !qp2)
2212                 goto error;
2213         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2214         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2215         cst1 = isl_upoly_as_cst(qp1->upoly);
2216         cst2 = isl_upoly_as_cst(qp2->upoly);
2217         cmp = isl_upoly_cmp(cst1, cst2);
2218
2219         if (cmp <= 0) {
2220                 isl_qpolynomial_free(qp2);
2221         } else {
2222                 isl_qpolynomial_free(qp1);
2223                 qp1 = qp2;
2224         }
2225         return qp1;
2226 error:
2227         isl_qpolynomial_free(qp1);
2228         isl_qpolynomial_free(qp2);
2229         return NULL;
2230 }
2231
2232 __isl_give isl_qpolynomial *isl_qpolynomial_max_cst(
2233         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2234 {
2235         struct isl_upoly_cst *cst1, *cst2;
2236         int cmp;
2237
2238         if (!qp1 || !qp2)
2239                 goto error;
2240         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2241         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2242         cst1 = isl_upoly_as_cst(qp1->upoly);
2243         cst2 = isl_upoly_as_cst(qp2->upoly);
2244         cmp = isl_upoly_cmp(cst1, cst2);
2245
2246         if (cmp >= 0) {
2247                 isl_qpolynomial_free(qp2);
2248         } else {
2249                 isl_qpolynomial_free(qp1);
2250                 qp1 = qp2;
2251         }
2252         return qp1;
2253 error:
2254         isl_qpolynomial_free(qp1);
2255         isl_qpolynomial_free(qp2);
2256         return NULL;
2257 }
2258
2259 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2260         __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2261         unsigned first, unsigned n)
2262 {
2263         unsigned total;
2264         unsigned g_pos;
2265         int *exp;
2266
2267         if (n == 0)
2268                 return qp;
2269
2270         qp = isl_qpolynomial_cow(qp);
2271         if (!qp)
2272                 return NULL;
2273
2274         isl_assert(qp->div->ctx, first <= isl_dim_size(qp->dim, type),
2275                     goto error);
2276
2277         g_pos = pos(qp->dim, type) + first;
2278
2279         qp->div = isl_mat_insert_cols(qp->div, 2 + g_pos, n);
2280         if (!qp->div)
2281                 goto error;
2282
2283         total = qp->div->n_col - 2;
2284         if (total > g_pos) {
2285                 int i;
2286                 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
2287                 if (!exp)
2288                         goto error;
2289                 for (i = 0; i < total - g_pos; ++i)
2290                         exp[i] = i + n;
2291                 qp->upoly = expand(qp->upoly, exp, g_pos);
2292                 free(exp);
2293                 if (!qp->upoly)
2294                         goto error;
2295         }
2296
2297         qp->dim = isl_dim_insert(qp->dim, type, first, n);
2298         if (!qp->dim)
2299                 goto error;
2300
2301         return qp;
2302 error:
2303         isl_qpolynomial_free(qp);
2304         return NULL;
2305 }
2306
2307 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
2308         __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
2309 {
2310         unsigned pos;
2311
2312         pos = isl_qpolynomial_dim(qp, type);
2313
2314         return isl_qpolynomial_insert_dims(qp, type, pos, n);
2315 }
2316
2317 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
2318         __isl_take isl_pw_qpolynomial *pwqp,
2319         enum isl_dim_type type, unsigned n)
2320 {
2321         unsigned pos;
2322
2323         pos = isl_pw_qpolynomial_dim(pwqp, type);
2324
2325         return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
2326 }
2327
2328 static int *reordering_move(isl_ctx *ctx,
2329         unsigned len, unsigned dst, unsigned src, unsigned n)
2330 {
2331         int i;
2332         int *reordering;
2333
2334         reordering = isl_alloc_array(ctx, int, len);
2335         if (!reordering)
2336                 return NULL;
2337
2338         if (dst <= src) {
2339                 for (i = 0; i < dst; ++i)
2340                         reordering[i] = i;
2341                 for (i = 0; i < n; ++i)
2342                         reordering[src + i] = dst + i;
2343                 for (i = 0; i < src - dst; ++i)
2344                         reordering[dst + i] = dst + n + i;
2345                 for (i = 0; i < len - src - n; ++i)
2346                         reordering[src + n + i] = src + n + i;
2347         } else {
2348                 for (i = 0; i < src; ++i)
2349                         reordering[i] = i;
2350                 for (i = 0; i < n; ++i)
2351                         reordering[src + i] = dst + i;
2352                 for (i = 0; i < dst - src; ++i)
2353                         reordering[src + n + i] = src + i;
2354                 for (i = 0; i < len - dst - n; ++i)
2355                         reordering[dst + n + i] = dst + n + i;
2356         }
2357
2358         return reordering;
2359 }
2360
2361 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
2362         __isl_take isl_qpolynomial *qp,
2363         enum isl_dim_type dst_type, unsigned dst_pos,
2364         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
2365 {
2366         unsigned g_dst_pos;
2367         unsigned g_src_pos;
2368         int *reordering;
2369
2370         qp = isl_qpolynomial_cow(qp);
2371         if (!qp)
2372                 return NULL;
2373
2374         isl_assert(qp->dim->ctx, src_pos + n <= isl_dim_size(qp->dim, src_type),
2375                 goto error);
2376
2377         g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
2378         g_src_pos = pos(qp->dim, src_type) + src_pos;
2379         if (dst_type > src_type)
2380                 g_dst_pos -= n;
2381
2382         qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
2383         if (!qp->div)
2384                 goto error;
2385         qp = sort_divs(qp);
2386         if (!qp)
2387                 goto error;
2388
2389         reordering = reordering_move(qp->dim->ctx,
2390                                 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
2391         if (!reordering)
2392                 goto error;
2393
2394         qp->upoly = reorder(qp->upoly, reordering);
2395         free(reordering);
2396         if (!qp->upoly)
2397                 goto error;
2398
2399         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
2400         if (!qp->dim)
2401                 goto error;
2402
2403         return qp;
2404 error:
2405         isl_qpolynomial_free(qp);
2406         return NULL;
2407 }
2408
2409 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_dim *dim,
2410         isl_int *f, isl_int denom)
2411 {
2412         struct isl_upoly *up;
2413
2414         if (!dim)
2415                 return NULL;
2416
2417         up = isl_upoly_from_affine(dim->ctx, f, denom, 1 + isl_dim_total(dim));
2418
2419         return isl_qpolynomial_alloc(dim, 0, up);
2420 }
2421
2422 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
2423         __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
2424 {
2425         isl_int denom;
2426         isl_dim *dim;
2427         struct isl_upoly *up;
2428         isl_qpolynomial *qp;
2429         int sgn;
2430
2431         if (!c)
2432                 return NULL;
2433
2434         isl_int_init(denom);
2435
2436         isl_constraint_get_coefficient(c, type, pos, &denom);
2437         isl_constraint_set_coefficient(c, type, pos, c->ctx->zero);
2438         sgn = isl_int_sgn(denom);
2439         isl_int_abs(denom, denom);
2440         up = isl_upoly_from_affine(c->ctx, c->line[0], denom,
2441                                         1 + isl_constraint_dim(c, isl_dim_all));
2442         if (sgn < 0)
2443                 isl_int_neg(denom, denom);
2444         isl_constraint_set_coefficient(c, type, pos, denom);
2445
2446         dim = isl_dim_copy(c->bmap->dim);
2447
2448         isl_int_clear(denom);
2449         isl_constraint_free(c);
2450
2451         qp = isl_qpolynomial_alloc(dim, 0, up);
2452         if (sgn > 0)
2453                 qp = isl_qpolynomial_neg(qp);
2454         return qp;
2455 }
2456
2457 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
2458  * in "qp" by subs[i].
2459  */
2460 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
2461         __isl_take isl_qpolynomial *qp,
2462         enum isl_dim_type type, unsigned first, unsigned n,
2463         __isl_keep isl_qpolynomial **subs)
2464 {
2465         int i;
2466         struct isl_upoly **ups;
2467
2468         if (n == 0)
2469                 return qp;
2470
2471         qp = isl_qpolynomial_cow(qp);
2472         if (!qp)
2473                 return NULL;
2474         for (i = 0; i < n; ++i)
2475                 if (!subs[i])
2476                         goto error;
2477
2478         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
2479                         goto error);
2480
2481         for (i = 0; i < n; ++i)
2482                 isl_assert(qp->dim->ctx, isl_dim_equal(qp->dim, subs[i]->dim),
2483                                 goto error);
2484
2485         isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
2486         for (i = 0; i < n; ++i)
2487                 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
2488
2489         first += pos(qp->dim, type);
2490
2491         ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
2492         if (!ups)
2493                 goto error;
2494         for (i = 0; i < n; ++i)
2495                 ups[i] = subs[i]->upoly;
2496
2497         qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
2498
2499         free(ups);
2500
2501         if (!qp->upoly)
2502                 goto error;
2503
2504         return qp;
2505 error:
2506         isl_qpolynomial_free(qp);
2507         return NULL;
2508 }
2509
2510 __isl_give isl_basic_set *add_div_constraints(__isl_take isl_basic_set *bset,
2511         __isl_take isl_mat *div)
2512 {
2513         int i;
2514         unsigned total;
2515
2516         if (!bset || !div)
2517                 goto error;
2518
2519         bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2520         if (!bset)
2521                 goto error;
2522         total = isl_basic_set_total_dim(bset);
2523         for (i = 0; i < div->n_row; ++i)
2524                 if (isl_basic_set_add_div_constraints_var(bset,
2525                                     total - div->n_row + i, div->row[i]) < 0)
2526                         goto error;
2527
2528         isl_mat_free(div);
2529         return bset;
2530 error:
2531         isl_mat_free(div);
2532         isl_basic_set_free(bset);
2533         return NULL;
2534 }
2535
2536 /* Extend "bset" with extra set dimensions for each integer division
2537  * in "qp" and then call "fn" with the extended bset and the polynomial
2538  * that results from replacing each of the integer divisions by the
2539  * corresponding extra set dimension.
2540  */
2541 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
2542         __isl_keep isl_basic_set *bset,
2543         int (*fn)(__isl_take isl_basic_set *bset,
2544                   __isl_take isl_qpolynomial *poly, void *user), void *user)
2545 {
2546         isl_dim *dim;
2547         isl_mat *div;
2548         isl_qpolynomial *poly;
2549
2550         if (!qp || !bset)
2551                 goto error;
2552         if (qp->div->n_row == 0)
2553                 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
2554                           user);
2555
2556         div = isl_mat_copy(qp->div);
2557         dim = isl_dim_copy(qp->dim);
2558         dim = isl_dim_add(dim, isl_dim_set, qp->div->n_row);
2559         poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
2560         bset = isl_basic_set_copy(bset);
2561         bset = isl_basic_set_add(bset, isl_dim_set, qp->div->n_row);
2562         bset = add_div_constraints(bset, div);
2563
2564         return fn(bset, poly, user);
2565 error:
2566         return -1;
2567 }
2568
2569 /* Return total degree in variables first (inclusive) up to last (exclusive).
2570  */
2571 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
2572 {
2573         int deg = -1;
2574         int i;
2575         struct isl_upoly_rec *rec;
2576
2577         if (!up)
2578                 return -2;
2579         if (isl_upoly_is_zero(up))
2580                 return -1;
2581         if (isl_upoly_is_cst(up) || up->var < first)
2582                 return 0;
2583
2584         rec = isl_upoly_as_rec(up);
2585         if (!rec)
2586                 return -2;
2587
2588         for (i = 0; i < rec->n; ++i) {
2589                 int d;
2590
2591                 if (isl_upoly_is_zero(rec->p[i]))
2592                         continue;
2593                 d = isl_upoly_degree(rec->p[i], first, last);
2594                 if (up->var < last)
2595                         d += i;
2596                 if (d > deg)
2597                         deg = d;
2598         }
2599
2600         return deg;
2601 }
2602
2603 /* Return total degree in set variables.
2604  */
2605 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
2606 {
2607         unsigned ovar;
2608         unsigned nvar;
2609
2610         if (!poly)
2611                 return -2;
2612
2613         ovar = isl_dim_offset(poly->dim, isl_dim_set);
2614         nvar = isl_dim_size(poly->dim, isl_dim_set);
2615         return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
2616 }
2617
2618 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
2619         unsigned pos, int deg)
2620 {
2621         int i;
2622         struct isl_upoly_rec *rec;
2623
2624         if (!up)
2625                 return NULL;
2626
2627         if (isl_upoly_is_cst(up) || up->var < pos) {
2628                 if (deg == 0)
2629                         return isl_upoly_copy(up);
2630                 else
2631                         return isl_upoly_zero(up->ctx);
2632         }
2633
2634         rec = isl_upoly_as_rec(up);
2635         if (!rec)
2636                 return NULL;
2637
2638         if (up->var == pos) {
2639                 if (deg < rec->n)
2640                         return isl_upoly_copy(rec->p[deg]);
2641                 else
2642                         return isl_upoly_zero(up->ctx);
2643         }
2644
2645         up = isl_upoly_copy(up);
2646         up = isl_upoly_cow(up);
2647         rec = isl_upoly_as_rec(up);
2648         if (!rec)
2649                 goto error;
2650
2651         for (i = 0; i < rec->n; ++i) {
2652                 struct isl_upoly *t;
2653                 t = isl_upoly_coeff(rec->p[i], pos, deg);
2654                 if (!t)
2655                         goto error;
2656                 isl_upoly_free(rec->p[i]);
2657                 rec->p[i] = t;
2658         }
2659
2660         return up;
2661 error:
2662         isl_upoly_free(up);
2663         return NULL;
2664 }
2665
2666 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
2667  */
2668 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
2669         __isl_keep isl_qpolynomial *qp,
2670         enum isl_dim_type type, unsigned t_pos, int deg)
2671 {
2672         unsigned g_pos;
2673         struct isl_upoly *up;
2674         isl_qpolynomial *c;
2675
2676         if (!qp)
2677                 return NULL;
2678
2679         isl_assert(qp->div->ctx, t_pos < isl_dim_size(qp->dim, type),
2680                         return NULL);
2681
2682         g_pos = pos(qp->dim, type) + t_pos;
2683         up = isl_upoly_coeff(qp->upoly, g_pos, deg);
2684
2685         c = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row, up);
2686         if (!c)
2687                 return NULL;
2688         isl_mat_free(c->div);
2689         c->div = isl_mat_copy(qp->div);
2690         if (!c->div)
2691                 goto error;
2692         return c;
2693 error:
2694         isl_qpolynomial_free(c);
2695         return NULL;
2696 }
2697
2698 /* Homogenize the polynomial in the variables first (inclusive) up to
2699  * last (exclusive) by inserting powers of variable first.
2700  * Variable first is assumed not to appear in the input.
2701  */
2702 __isl_give struct isl_upoly *isl_upoly_homogenize(
2703         __isl_take struct isl_upoly *up, int deg, int target,
2704         int first, int last)
2705 {
2706         int i;
2707         struct isl_upoly_rec *rec;
2708
2709         if (!up)
2710                 return NULL;
2711         if (isl_upoly_is_zero(up))
2712                 return up;
2713         if (deg == target)
2714                 return up;
2715         if (isl_upoly_is_cst(up) || up->var < first) {
2716                 struct isl_upoly *hom;
2717
2718                 hom = isl_upoly_pow(up->ctx, first, target - deg);
2719                 if (!hom)
2720                         goto error;
2721                 rec = isl_upoly_as_rec(hom);
2722                 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
2723
2724                 return hom;
2725         }
2726
2727         up = isl_upoly_cow(up);
2728         rec = isl_upoly_as_rec(up);
2729         if (!rec)
2730                 goto error;
2731
2732         for (i = 0; i < rec->n; ++i) {
2733                 if (isl_upoly_is_zero(rec->p[i]))
2734                         continue;
2735                 rec->p[i] = isl_upoly_homogenize(rec->p[i],
2736                                 up->var < last ? deg + i : i, target,
2737                                 first, last);
2738                 if (!rec->p[i])
2739                         goto error;
2740         }
2741
2742         return up;
2743 error:
2744         isl_upoly_free(up);
2745         return NULL;
2746 }
2747
2748 /* Homogenize the polynomial in the set variables by introducing
2749  * powers of an extra set variable at position 0.
2750  */
2751 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
2752         __isl_take isl_qpolynomial *poly)
2753 {
2754         unsigned ovar;
2755         unsigned nvar;
2756         int deg = isl_qpolynomial_degree(poly);
2757
2758         if (deg < -1)
2759                 goto error;
2760
2761         poly = isl_qpolynomial_insert_dims(poly, isl_dim_set, 0, 1);
2762         poly = isl_qpolynomial_cow(poly);
2763         if (!poly)
2764                 goto error;
2765
2766         ovar = isl_dim_offset(poly->dim, isl_dim_set);
2767         nvar = isl_dim_size(poly->dim, isl_dim_set);
2768         poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
2769                                                 ovar, ovar + nvar);
2770         if (!poly->upoly)
2771                 goto error;
2772
2773         return poly;
2774 error:
2775         isl_qpolynomial_free(poly);
2776         return NULL;
2777 }
2778
2779 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
2780         __isl_take isl_mat *div)
2781 {
2782         isl_term *term;
2783         int n;
2784
2785         if (!dim || !div)
2786                 goto error;
2787
2788         n = isl_dim_total(dim) + div->n_row;
2789
2790         term = isl_calloc(dim->ctx, struct isl_term,
2791                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
2792         if (!term)
2793                 goto error;
2794
2795         term->ref = 1;
2796         term->dim = dim;
2797         term->div = div;
2798         isl_int_init(term->n);
2799         isl_int_init(term->d);
2800
2801         return term;
2802 error:
2803         isl_dim_free(dim);
2804         isl_mat_free(div);
2805         return NULL;
2806 }
2807
2808 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
2809 {
2810         if (!term)
2811                 return NULL;
2812
2813         term->ref++;
2814         return term;
2815 }
2816
2817 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
2818 {
2819         int i;
2820         isl_term *dup;
2821         unsigned total;
2822
2823         if (term)
2824                 return NULL;
2825
2826         total = isl_dim_total(term->dim) + term->div->n_row;
2827
2828         dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
2829         if (!dup)
2830                 return NULL;
2831
2832         isl_int_set(dup->n, term->n);
2833         isl_int_set(dup->d, term->d);
2834
2835         for (i = 0; i < total; ++i)
2836                 dup->pow[i] = term->pow[i];
2837
2838         return dup;
2839 }
2840
2841 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
2842 {
2843         if (!term)
2844                 return NULL;
2845
2846         if (term->ref == 1)
2847                 return term;
2848         term->ref--;
2849         return isl_term_dup(term);
2850 }
2851
2852 void isl_term_free(__isl_take isl_term *term)
2853 {
2854         if (!term)
2855                 return;
2856
2857         if (--term->ref > 0)
2858                 return;
2859
2860         isl_dim_free(term->dim);
2861         isl_mat_free(term->div);
2862         isl_int_clear(term->n);
2863         isl_int_clear(term->d);
2864         free(term);
2865 }
2866
2867 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
2868 {
2869         if (!term)
2870                 return 0;
2871
2872         switch (type) {
2873         case isl_dim_param:
2874         case isl_dim_in:
2875         case isl_dim_out:       return isl_dim_size(term->dim, type);
2876         case isl_dim_div:       return term->div->n_row;
2877         case isl_dim_all:       return isl_dim_total(term->dim) + term->div->n_row;
2878         default:                return 0;
2879         }
2880 }
2881
2882 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
2883 {
2884         return term ? term->dim->ctx : NULL;
2885 }
2886
2887 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
2888 {
2889         if (!term)
2890                 return;
2891         isl_int_set(*n, term->n);
2892 }
2893
2894 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
2895 {
2896         if (!term)
2897                 return;
2898         isl_int_set(*d, term->d);
2899 }
2900
2901 int isl_term_get_exp(__isl_keep isl_term *term,
2902         enum isl_dim_type type, unsigned pos)
2903 {
2904         if (!term)
2905                 return -1;
2906
2907         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
2908
2909         if (type >= isl_dim_set)
2910                 pos += isl_dim_size(term->dim, isl_dim_param);
2911         if (type >= isl_dim_div)
2912                 pos += isl_dim_size(term->dim, isl_dim_set);
2913
2914         return term->pow[pos];
2915 }
2916
2917 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
2918 {
2919         isl_basic_map *bmap;
2920         unsigned total;
2921         int k;
2922
2923         if (!term)
2924                 return NULL;
2925
2926         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
2927                         return NULL);
2928
2929         total = term->div->n_col - term->div->n_row - 2;
2930         /* No nested divs for now */
2931         isl_assert(term->dim->ctx,
2932                 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
2933                                         term->div->n_row) == -1,
2934                 return NULL);
2935
2936         bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
2937         if ((k = isl_basic_map_alloc_div(bmap)) < 0)
2938                 goto error;
2939
2940         isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
2941
2942         return isl_basic_map_div(bmap, k);
2943 error:
2944         isl_basic_map_free(bmap);
2945         return NULL;
2946 }
2947
2948 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
2949         int (*fn)(__isl_take isl_term *term, void *user),
2950         __isl_take isl_term *term, void *user)
2951 {
2952         int i;
2953         struct isl_upoly_rec *rec;
2954
2955         if (!up || !term)
2956                 goto error;
2957
2958         if (isl_upoly_is_zero(up))
2959                 return term;
2960
2961         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
2962         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
2963         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
2964
2965         if (isl_upoly_is_cst(up)) {
2966                 struct isl_upoly_cst *cst;
2967                 cst = isl_upoly_as_cst(up);
2968                 if (!cst)
2969                         goto error;
2970                 term = isl_term_cow(term);
2971                 if (!term)
2972                         goto error;
2973                 isl_int_set(term->n, cst->n);
2974                 isl_int_set(term->d, cst->d);
2975                 if (fn(isl_term_copy(term), user) < 0)
2976                         goto error;
2977                 return term;
2978         }
2979
2980         rec = isl_upoly_as_rec(up);
2981         if (!rec)
2982                 goto error;
2983
2984         for (i = 0; i < rec->n; ++i) {
2985                 term = isl_term_cow(term);
2986                 if (!term)
2987                         goto error;
2988                 term->pow[up->var] = i;
2989                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2990                 if (!term)
2991                         goto error;
2992         }
2993         term->pow[up->var] = 0;
2994
2995         return term;
2996 error:
2997         isl_term_free(term);
2998         return NULL;
2999 }
3000
3001 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3002         int (*fn)(__isl_take isl_term *term, void *user), void *user)
3003 {
3004         isl_term *term;
3005
3006         if (!qp)
3007                 return -1;
3008
3009         term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
3010         if (!term)
3011                 return -1;
3012
3013         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3014
3015         isl_term_free(term);
3016
3017         return term ? 0 : -1;
3018 }
3019
3020 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3021 {
3022         struct isl_upoly *up;
3023         isl_qpolynomial *qp;
3024         int i, n;
3025
3026         if (!term)
3027                 return NULL;
3028
3029         n = isl_dim_total(term->dim) + term->div->n_row;
3030
3031         up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3032         for (i = 0; i < n; ++i) {
3033                 if (!term->pow[i])
3034                         continue;
3035                 up = isl_upoly_mul(up,
3036                         isl_upoly_pow(term->dim->ctx, i, term->pow[i]));
3037         }
3038
3039         qp = isl_qpolynomial_alloc(isl_dim_copy(term->dim), term->div->n_row, up);
3040         if (!qp)
3041                 goto error;
3042         isl_mat_free(qp->div);
3043         qp->div = isl_mat_copy(term->div);
3044         if (!qp->div)
3045                 goto error;
3046
3047         isl_term_free(term);
3048         return qp;
3049 error:
3050         isl_qpolynomial_free(qp);
3051         isl_term_free(term);
3052         return NULL;
3053 }
3054
3055 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3056         __isl_take isl_dim *dim)
3057 {
3058         int i;
3059         int extra;
3060         unsigned total;
3061
3062         if (!qp || !dim)
3063                 goto error;
3064
3065         if (isl_dim_equal(qp->dim, dim)) {
3066                 isl_dim_free(dim);
3067                 return qp;
3068         }
3069
3070         qp = isl_qpolynomial_cow(qp);
3071         if (!qp)
3072                 goto error;
3073
3074         extra = isl_dim_size(dim, isl_dim_set) -
3075                         isl_dim_size(qp->dim, isl_dim_set);
3076         total = isl_dim_total(qp->dim);
3077         if (qp->div->n_row) {
3078                 int *exp;
3079
3080                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3081                 if (!exp)
3082                         goto error;
3083                 for (i = 0; i < qp->div->n_row; ++i)
3084                         exp[i] = extra + i;
3085                 qp->upoly = expand(qp->upoly, exp, total);
3086                 free(exp);
3087                 if (!qp->upoly)
3088                         goto error;
3089         }
3090         qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3091         if (!qp->div)
3092                 goto error;
3093         for (i = 0; i < qp->div->n_row; ++i)
3094                 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3095
3096         isl_dim_free(qp->dim);
3097         qp->dim = dim;
3098
3099         return qp;
3100 error:
3101         isl_dim_free(dim);
3102         isl_qpolynomial_free(qp);
3103         return NULL;
3104 }
3105
3106 /* For each parameter or variable that does not appear in qp,
3107  * first eliminate the variable from all constraints and then set it to zero.
3108  */
3109 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3110         __isl_keep isl_qpolynomial *qp)
3111 {
3112         int *active = NULL;
3113         int i;
3114         int d;
3115         unsigned nparam;
3116         unsigned nvar;
3117
3118         if (!set || !qp)
3119                 goto error;
3120
3121         d = isl_dim_total(set->dim);
3122         active = isl_calloc_array(set->ctx, int, d);
3123         if (set_active(qp, active) < 0)
3124                 goto error;
3125
3126         for (i = 0; i < d; ++i)
3127                 if (!active[i])
3128                         break;
3129
3130         if (i == d) {
3131                 free(active);
3132                 return set;
3133         }
3134
3135         nparam = isl_dim_size(set->dim, isl_dim_param);
3136         nvar = isl_dim_size(set->dim, isl_dim_set);
3137         for (i = 0; i < nparam; ++i) {
3138                 if (active[i])
3139                         continue;
3140                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3141                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3142         }
3143         for (i = 0; i < nvar; ++i) {
3144                 if (active[nparam + i])
3145                         continue;
3146                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3147                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3148         }
3149
3150         free(active);
3151
3152         return set;
3153 error:
3154         free(active);
3155         isl_set_free(set);
3156         return NULL;
3157 }
3158
3159 struct isl_opt_data {
3160         isl_qpolynomial *qp;
3161         int first;
3162         isl_qpolynomial *opt;
3163         int max;
3164 };
3165
3166 static int opt_fn(__isl_take isl_point *pnt, void *user)
3167 {
3168         struct isl_opt_data *data = (struct isl_opt_data *)user;
3169         isl_qpolynomial *val;
3170
3171         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3172         if (data->first) {
3173                 data->first = 0;
3174                 data->opt = val;
3175         } else if (data->max) {
3176                 data->opt = isl_qpolynomial_max_cst(data->opt, val);
3177         } else {
3178                 data->opt = isl_qpolynomial_min_cst(data->opt, val);
3179         }
3180
3181         return 0;
3182 }
3183
3184 __isl_give isl_qpolynomial *isl_qpolynomial_opt_on_domain(
3185         __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3186 {
3187         struct isl_opt_data data = { NULL, 1, NULL, max };
3188
3189         if (!set || !qp)
3190                 goto error;
3191
3192         if (isl_upoly_is_cst(qp->upoly)) {
3193                 isl_set_free(set);
3194                 return qp;
3195         }
3196
3197         set = fix_inactive(set, qp);
3198
3199         data.qp = qp;
3200         if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3201                 goto error;
3202
3203         if (data.first)
3204                 data.opt = isl_qpolynomial_zero(isl_qpolynomial_get_dim(qp));
3205
3206         isl_set_free(set);
3207         isl_qpolynomial_free(qp);
3208         return data.opt;
3209 error:
3210         isl_set_free(set);
3211         isl_qpolynomial_free(qp);
3212         isl_qpolynomial_free(data.opt);
3213         return NULL;
3214 }
3215
3216 __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp,
3217         __isl_take isl_morph *morph)
3218 {
3219         int i;
3220         isl_ctx *ctx;
3221         struct isl_upoly *up;
3222         unsigned n_div;
3223         struct isl_upoly **subs;
3224         isl_mat *mat;
3225
3226         qp = isl_qpolynomial_cow(qp);
3227         if (!qp || !morph)
3228                 goto error;
3229
3230         ctx = qp->dim->ctx;
3231         isl_assert(ctx, isl_dim_equal(qp->dim, morph->dom->dim), goto error);
3232
3233         subs = isl_calloc_array(ctx, struct isl_upoly *, morph->inv->n_row - 1);
3234         if (!subs)
3235                 goto error;
3236
3237         for (i = 0; 1 + i < morph->inv->n_row; ++i)
3238                 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
3239                                         morph->inv->row[0][0], morph->inv->n_col);
3240
3241         qp->upoly = isl_upoly_subs(qp->upoly, 0, morph->inv->n_row - 1, subs);
3242
3243         for (i = 0; 1 + i < morph->inv->n_row; ++i)
3244                 isl_upoly_free(subs[i]);
3245         free(subs);
3246
3247         mat = isl_mat_diagonal(isl_mat_identity(ctx, 1), isl_mat_copy(morph->inv));
3248         mat = isl_mat_diagonal(mat, isl_mat_identity(ctx, qp->div->n_row));
3249         qp->div = isl_mat_product(qp->div, mat);
3250         isl_dim_free(qp->dim);
3251         qp->dim = isl_dim_copy(morph->ran->dim);
3252
3253         if (!qp->upoly || !qp->div || !qp->dim)
3254                 goto error;
3255
3256         isl_morph_free(morph);
3257
3258         return qp;
3259 error:
3260         isl_qpolynomial_free(qp);
3261         isl_morph_free(morph);
3262         return NULL;
3263 }
3264
3265 static int neg_entry(void **entry, void *user)
3266 {
3267         isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
3268
3269         *pwqp = isl_pw_qpolynomial_neg(*pwqp);
3270
3271         return *pwqp ? 0 : -1;
3272 }
3273
3274 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
3275         __isl_take isl_union_pw_qpolynomial *upwqp)
3276 {
3277         upwqp = isl_union_pw_qpolynomial_cow(upwqp);
3278         if (!upwqp)
3279                 return NULL;
3280
3281         if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
3282                                    &neg_entry, NULL) < 0)
3283                 goto error;
3284
3285         return upwqp;
3286 error:
3287         isl_union_pw_qpolynomial_free(upwqp);
3288         return NULL;
3289 }
3290
3291 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sub(
3292         __isl_take isl_union_pw_qpolynomial *upwqp1,
3293         __isl_take isl_union_pw_qpolynomial *upwqp2)
3294 {
3295         return isl_union_pw_qpolynomial_add(upwqp1,
3296                                         isl_union_pw_qpolynomial_neg(upwqp2));
3297 }
3298
3299 static int mul_entry(void **entry, void *user)
3300 {
3301         struct isl_union_pw_qpolynomial_match_bin_data *data = user;
3302         uint32_t hash;
3303         struct isl_hash_table_entry *entry2;
3304         isl_pw_qpolynomial *pwpq = *entry;
3305         int empty;
3306
3307         hash = isl_dim_get_hash(pwpq->dim);
3308         entry2 = isl_hash_table_find(data->u2->dim->ctx, &data->u2->table,
3309                                      hash, &has_dim, pwpq->dim, 0);
3310         if (!entry2)
3311                 return 0;
3312
3313         pwpq = isl_pw_qpolynomial_copy(pwpq);
3314         pwpq = isl_pw_qpolynomial_mul(pwpq,
3315                                       isl_pw_qpolynomial_copy(entry2->data));
3316
3317         empty = isl_pw_qpolynomial_is_zero(pwpq);
3318         if (empty < 0) {
3319                 isl_pw_qpolynomial_free(pwpq);
3320                 return -1;
3321         }
3322         if (empty) {
3323                 isl_pw_qpolynomial_free(pwpq);
3324                 return 0;
3325         }
3326
3327         data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(data->res, pwpq);
3328
3329         return 0;
3330 }
3331
3332 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
3333         __isl_take isl_union_pw_qpolynomial *upwqp1,
3334         __isl_take isl_union_pw_qpolynomial *upwqp2)
3335 {
3336         return match_bin_op(upwqp1, upwqp2, &mul_entry);
3337 }
3338
3339 /* Reorder the columns of the given div definitions according to the
3340  * given reordering.
3341  */
3342 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
3343         __isl_take isl_reordering *r)
3344 {
3345         int i, j;
3346         isl_mat *mat;
3347         int extra;
3348
3349         if (!div || !r)
3350                 goto error;
3351
3352         extra = isl_dim_total(r->dim) + div->n_row - r->len;
3353         mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
3354         if (!mat)
3355                 goto error;
3356
3357         for (i = 0; i < div->n_row; ++i) {
3358                 isl_seq_cpy(mat->row[i], div->row[i], 2);
3359                 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
3360                 for (j = 0; j < r->len; ++j)
3361                         isl_int_set(mat->row[i][2 + r->pos[j]],
3362                                     div->row[i][2 + j]);
3363         }
3364
3365         isl_reordering_free(r);
3366         isl_mat_free(div);
3367         return mat;
3368 error:
3369         isl_reordering_free(r);
3370         isl_mat_free(div);
3371         return NULL;
3372 }
3373
3374 /* Reorder the dimension of "qp" according to the given reordering.
3375  */
3376 __isl_give isl_qpolynomial *isl_qpolynomial_realign(
3377         __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
3378 {
3379         qp = isl_qpolynomial_cow(qp);
3380         if (!qp)
3381                 goto error;
3382
3383         r = isl_reordering_extend(r, qp->div->n_row);
3384         if (!r)
3385                 goto error;
3386
3387         qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
3388         if (!qp->div)
3389                 goto error;
3390
3391         qp->upoly = reorder(qp->upoly, r->pos);
3392         if (!qp->upoly)
3393                 goto error;
3394
3395         qp = isl_qpolynomial_reset_dim(qp, isl_dim_copy(r->dim));
3396
3397         isl_reordering_free(r);
3398         return qp;
3399 error:
3400         isl_qpolynomial_free(qp);
3401         isl_reordering_free(r);
3402         return NULL;
3403 }