isl_polynomial.c

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