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