isl_polynomial.c

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