isl_polynomial.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the MIT 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 #define ISL_DIM_H
  13 #include <isl_ctx_private.h>
  14 #include <isl_map_private.h>
  15 #include <isl_factorization.h>
  16 #include <isl/lp.h>
  17 #include <isl/seq.h>
  18 #include <isl_union_map_private.h>
  19 #include <isl_constraint_private.h>
  20 #include <isl_polynomial_private.h>
  21 #include <isl_point_private.h>
  22 #include <isl_space_private.h>
  23 #include <isl_mat_private.h>
  24 #include <isl_range.h>
  25 #include <isl_local_space_private.h>
  26 #include <isl_aff_private.h>
  27 #include <isl_val_private.h>
  28 #include <isl_config.h>
  29
  30 static unsigned pos(__isl_keep isl_space *dim, enum isl_dim_type type)
  31 {
  32         switch (type) {
  33         case isl_dim_param:     return 0;
  34         case isl_dim_in:        return dim->nparam;
  35         case isl_dim_out:       return dim->nparam + dim->n_in;
  36         default:                return 0;
  37         }
  38 }
  39
  40 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
  41 {
  42         if (!up)
  43                 return -1;
  44
  45         return up->var < 0;
  46 }
  47
  48 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
  49 {
  50         if (!up)
  51                 return NULL;
  52
  53         isl_assert(up->ctx, up->var < 0, return NULL);
  54
  55         return (struct isl_upoly_cst *)up;
  56 }
  57
  58 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
  59 {
  60         if (!up)
  61                 return NULL;
  62
  63         isl_assert(up->ctx, up->var >= 0, return NULL);
  64
  65         return (struct isl_upoly_rec *)up;
  66 }
  67
  68 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
  69         __isl_keep struct isl_upoly *up2)
  70 {
  71         int i;
  72         struct isl_upoly_rec *rec1, *rec2;
  73
  74         if (!up1 || !up2)
  75                 return -1;
  76         if (up1 == up2)
  77                 return 1;
  78         if (up1->var != up2->var)
  79                 return 0;
  80         if (isl_upoly_is_cst(up1)) {
  81                 struct isl_upoly_cst *cst1, *cst2;
  82                 cst1 = isl_upoly_as_cst(up1);
  83                 cst2 = isl_upoly_as_cst(up2);
  84                 if (!cst1 || !cst2)
  85                         return -1;
  86                 return isl_int_eq(cst1->n, cst2->n) &&
  87                        isl_int_eq(cst1->d, cst2->d);
  88         }
  89
  90         rec1 = isl_upoly_as_rec(up1);
  91         rec2 = isl_upoly_as_rec(up2);
  92         if (!rec1 || !rec2)
  93                 return -1;
  94
  95         if (rec1->n != rec2->n)
  96                 return 0;
  97
  98         for (i = 0; i < rec1->n; ++i) {
  99                 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
 100                 if (eq < 0 || !eq)
 101                         return eq;
 102         }
 103
 104         return 1;
 105 }
 106
 107 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
 108 {
 109         struct isl_upoly_cst *cst;
 110
 111         if (!up)
 112                 return -1;
 113         if (!isl_upoly_is_cst(up))
 114                 return 0;
 115
 116         cst = isl_upoly_as_cst(up);
 117         if (!cst)
 118                 return -1;
 119
 120         return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
 121 }
 122
 123 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
 124 {
 125         struct isl_upoly_cst *cst;
 126
 127         if (!up)
 128                 return 0;
 129         if (!isl_upoly_is_cst(up))
 130                 return 0;
 131
 132         cst = isl_upoly_as_cst(up);
 133         if (!cst)
 134                 return 0;
 135
 136         return isl_int_sgn(cst->n);
 137 }
 138
 139 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
 140 {
 141         struct isl_upoly_cst *cst;
 142
 143         if (!up)
 144                 return -1;
 145         if (!isl_upoly_is_cst(up))
 146                 return 0;
 147
 148         cst = isl_upoly_as_cst(up);
 149         if (!cst)
 150                 return -1;
 151
 152         return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
 153 }
 154
 155 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
 156 {
 157         struct isl_upoly_cst *cst;
 158
 159         if (!up)
 160                 return -1;
 161         if (!isl_upoly_is_cst(up))
 162                 return 0;
 163
 164         cst = isl_upoly_as_cst(up);
 165         if (!cst)
 166                 return -1;
 167
 168         return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
 169 }
 170
 171 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
 172 {
 173         struct isl_upoly_cst *cst;
 174
 175         if (!up)
 176                 return -1;
 177         if (!isl_upoly_is_cst(up))
 178                 return 0;
 179
 180         cst = isl_upoly_as_cst(up);
 181         if (!cst)
 182                 return -1;
 183
 184         return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
 185 }
 186
 187 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
 188 {
 189         struct isl_upoly_cst *cst;
 190
 191         if (!up)
 192                 return -1;
 193         if (!isl_upoly_is_cst(up))
 194                 return 0;
 195
 196         cst = isl_upoly_as_cst(up);
 197         if (!cst)
 198                 return -1;
 199
 200         return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
 201 }
 202
 203 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
 204 {
 205         struct isl_upoly_cst *cst;
 206
 207         if (!up)
 208                 return -1;
 209         if (!isl_upoly_is_cst(up))
 210                 return 0;
 211
 212         cst = isl_upoly_as_cst(up);
 213         if (!cst)
 214                 return -1;
 215
 216         return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
 217 }
 218
 219 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
 220 {
 221         struct isl_upoly_cst *cst;
 222
 223         cst = isl_alloc_type(ctx, struct isl_upoly_cst);
 224         if (!cst)
 225                 return NULL;
 226
 227         cst->up.ref = 1;
 228         cst->up.ctx = ctx;
 229         isl_ctx_ref(ctx);
 230         cst->up.var = -1;
 231
 232         isl_int_init(cst->n);
 233         isl_int_init(cst->d);
 234
 235         return cst;
 236 }
 237
 238 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
 239 {
 240         struct isl_upoly_cst *cst;
 241
 242         cst = isl_upoly_cst_alloc(ctx);
 243         if (!cst)
 244                 return NULL;
 245
 246         isl_int_set_si(cst->n, 0);
 247         isl_int_set_si(cst->d, 1);
 248
 249         return &cst->up;
 250 }
 251
 252 __isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
 253 {
 254         struct isl_upoly_cst *cst;
 255
 256         cst = isl_upoly_cst_alloc(ctx);
 257         if (!cst)
 258                 return NULL;
 259
 260         isl_int_set_si(cst->n, 1);
 261         isl_int_set_si(cst->d, 1);
 262
 263         return &cst->up;
 264 }
 265
 266 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
 267 {
 268         struct isl_upoly_cst *cst;
 269
 270         cst = isl_upoly_cst_alloc(ctx);
 271         if (!cst)
 272                 return NULL;
 273
 274         isl_int_set_si(cst->n, 1);
 275         isl_int_set_si(cst->d, 0);
 276
 277         return &cst->up;
 278 }
 279
 280 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
 281 {
 282         struct isl_upoly_cst *cst;
 283
 284         cst = isl_upoly_cst_alloc(ctx);
 285         if (!cst)
 286                 return NULL;
 287
 288         isl_int_set_si(cst->n, -1);
 289         isl_int_set_si(cst->d, 0);
 290
 291         return &cst->up;
 292 }
 293
 294 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
 295 {
 296         struct isl_upoly_cst *cst;
 297
 298         cst = isl_upoly_cst_alloc(ctx);
 299         if (!cst)
 300                 return NULL;
 301
 302         isl_int_set_si(cst->n, 0);
 303         isl_int_set_si(cst->d, 0);
 304
 305         return &cst->up;
 306 }
 307
 308 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
 309         isl_int n, isl_int d)
 310 {
 311         struct isl_upoly_cst *cst;
 312
 313         cst = isl_upoly_cst_alloc(ctx);
 314         if (!cst)
 315                 return NULL;
 316
 317         isl_int_set(cst->n, n);
 318         isl_int_set(cst->d, d);
 319
 320         return &cst->up;
 321 }
 322
 323 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
 324         int var, int size)
 325 {
 326         struct isl_upoly_rec *rec;
 327
 328         isl_assert(ctx, var >= 0, return NULL);
 329         isl_assert(ctx, size >= 0, return NULL);
 330         rec = isl_calloc(ctx, struct isl_upoly_rec,
 331                         sizeof(struct isl_upoly_rec) +
 332                         size * sizeof(struct isl_upoly *));
 333         if (!rec)
 334                 return NULL;
 335
 336         rec->up.ref = 1;
 337         rec->up.ctx = ctx;
 338         isl_ctx_ref(ctx);
 339         rec->up.var = var;
 340
 341         rec->n = 0;
 342         rec->size = size;
 343
 344         return rec;
 345 }
 346
 347 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
 348         __isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
 349 {
 350         qp = isl_qpolynomial_cow(qp);
 351         if (!qp || !dim)
 352                 goto error;
 353
 354         isl_space_free(qp->dim);
 355         qp->dim = dim;
 356
 357         return qp;
 358 error:
 359         isl_qpolynomial_free(qp);
 360         isl_space_free(dim);
 361         return NULL;
 362 }
 363
 364 /* Reset the space of "qp".  This function is called from isl_pw_templ.c
 365  * and doesn't know if the space of an element object is represented
 366  * directly or through its domain.  It therefore passes along both.
 367  */
 368 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
 369         __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
 370         __isl_take isl_space *domain)
 371 {
 372         isl_space_free(space);
 373         return isl_qpolynomial_reset_domain_space(qp, domain);
 374 }
 375
 376 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
 377 {
 378         return qp ? qp->dim->ctx : NULL;
 379 }
 380
 381 __isl_give isl_space *isl_qpolynomial_get_domain_space(
 382         __isl_keep isl_qpolynomial *qp)
 383 {
 384         return qp ? isl_space_copy(qp->dim) : NULL;
 385 }
 386
 387 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
 388 {
 389         isl_space *space;
 390         if (!qp)
 391                 return NULL;
 392         space = isl_space_copy(qp->dim);
 393         space = isl_space_from_domain(space);
 394         space = isl_space_add_dims(space, isl_dim_out, 1);
 395         return space;
 396 }
 397
 398 /* Externally, an isl_qpolynomial has a map space, but internally, the
 399  * ls field corresponds to the domain of that space.
 400  */
 401 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
 402         enum isl_dim_type type)
 403 {
 404         if (!qp)
 405                 return 0;
 406         if (type == isl_dim_out)
 407                 return 1;
 408         if (type == isl_dim_in)
 409                 type = isl_dim_set;
 410         return isl_space_dim(qp->dim, type);
 411 }
 412
 413 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
 414 {
 415         return qp ? isl_upoly_is_zero(qp->upoly) : -1;
 416 }
 417
 418 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
 419 {
 420         return qp ? isl_upoly_is_one(qp->upoly) : -1;
 421 }
 422
 423 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
 424 {
 425         return qp ? isl_upoly_is_nan(qp->upoly) : -1;
 426 }
 427
 428 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
 429 {
 430         return qp ? isl_upoly_is_infty(qp->upoly) : -1;
 431 }
 432
 433 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
 434 {
 435         return qp ? isl_upoly_is_neginfty(qp->upoly) : -1;
 436 }
 437
 438 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
 439 {
 440         return qp ? isl_upoly_sgn(qp->upoly) : 0;
 441 }
 442
 443 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
 444 {
 445         isl_int_clear(cst->n);
 446         isl_int_clear(cst->d);
 447 }
 448
 449 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
 450 {
 451         int i;
 452
 453         for (i = 0; i < rec->n; ++i)
 454                 isl_upoly_free(rec->p[i]);
 455 }
 456
 457 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
 458 {
 459         if (!up)
 460                 return NULL;
 461
 462         up->ref++;
 463         return up;
 464 }
 465
 466 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
 467 {
 468         struct isl_upoly_cst *cst;
 469         struct isl_upoly_cst *dup;
 470
 471         cst = isl_upoly_as_cst(up);
 472         if (!cst)
 473                 return NULL;
 474
 475         dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
 476         if (!dup)
 477                 return NULL;
 478         isl_int_set(dup->n, cst->n);
 479         isl_int_set(dup->d, cst->d);
 480
 481         return &dup->up;
 482 }
 483
 484 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
 485 {
 486         int i;
 487         struct isl_upoly_rec *rec;
 488         struct isl_upoly_rec *dup;
 489
 490         rec = isl_upoly_as_rec(up);
 491         if (!rec)
 492                 return NULL;
 493
 494         dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
 495         if (!dup)
 496                 return NULL;
 497
 498         for (i = 0; i < rec->n; ++i) {
 499                 dup->p[i] = isl_upoly_copy(rec->p[i]);
 500                 if (!dup->p[i])
 501                         goto error;
 502                 dup->n++;
 503         }
 504
 505         return &dup->up;
 506 error:
 507         isl_upoly_free(&dup->up);
 508         return NULL;
 509 }
 510
 511 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
 512 {
 513         if (!up)
 514                 return NULL;
 515
 516         if (isl_upoly_is_cst(up))
 517                 return isl_upoly_dup_cst(up);
 518         else
 519                 return isl_upoly_dup_rec(up);
 520 }
 521
 522 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
 523 {
 524         if (!up)
 525                 return NULL;
 526
 527         if (up->ref == 1)
 528                 return up;
 529         up->ref--;
 530         return isl_upoly_dup(up);
 531 }
 532
 533 void isl_upoly_free(__isl_take struct isl_upoly *up)
 534 {
 535         if (!up)
 536                 return;
 537
 538         if (--up->ref > 0)
 539                 return;
 540
 541         if (up->var < 0)
 542                 upoly_free_cst((struct isl_upoly_cst *)up);
 543         else
 544                 upoly_free_rec((struct isl_upoly_rec *)up);
 545
 546         isl_ctx_deref(up->ctx);
 547         free(up);
 548 }
 549
 550 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
 551 {
 552         isl_int gcd;
 553
 554         isl_int_init(gcd);
 555         isl_int_gcd(gcd, cst->n, cst->d);
 556         if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
 557                 isl_int_divexact(cst->n, cst->n, gcd);
 558                 isl_int_divexact(cst->d, cst->d, gcd);
 559         }
 560         isl_int_clear(gcd);
 561 }
 562
 563 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
 564         __isl_take struct isl_upoly *up2)
 565 {
 566         struct isl_upoly_cst *cst1;
 567         struct isl_upoly_cst *cst2;
 568
 569         up1 = isl_upoly_cow(up1);
 570         if (!up1 || !up2)
 571                 goto error;
 572
 573         cst1 = isl_upoly_as_cst(up1);
 574         cst2 = isl_upoly_as_cst(up2);
 575
 576         if (isl_int_eq(cst1->d, cst2->d))
 577                 isl_int_add(cst1->n, cst1->n, cst2->n);
 578         else {
 579                 isl_int_mul(cst1->n, cst1->n, cst2->d);
 580                 isl_int_addmul(cst1->n, cst2->n, cst1->d);
 581                 isl_int_mul(cst1->d, cst1->d, cst2->d);
 582         }
 583
 584         isl_upoly_cst_reduce(cst1);
 585
 586         isl_upoly_free(up2);
 587         return up1;
 588 error:
 589         isl_upoly_free(up1);
 590         isl_upoly_free(up2);
 591         return NULL;
 592 }
 593
 594 static __isl_give struct isl_upoly *replace_by_zero(
 595         __isl_take struct isl_upoly *up)
 596 {
 597         struct isl_ctx *ctx;
 598
 599         if (!up)
 600                 return NULL;
 601         ctx = up->ctx;
 602         isl_upoly_free(up);
 603         return isl_upoly_zero(ctx);
 604 }
 605
 606 static __isl_give struct isl_upoly *replace_by_constant_term(
 607         __isl_take struct isl_upoly *up)
 608 {
 609         struct isl_upoly_rec *rec;
 610         struct isl_upoly *cst;
 611
 612         if (!up)
 613                 return NULL;
 614
 615         rec = isl_upoly_as_rec(up);
 616         if (!rec)
 617                 goto error;
 618         cst = isl_upoly_copy(rec->p[0]);
 619         isl_upoly_free(up);
 620         return cst;
 621 error:
 622         isl_upoly_free(up);
 623         return NULL;
 624 }
 625
 626 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
 627         __isl_take struct isl_upoly *up2)
 628 {
 629         int i;
 630         struct isl_upoly_rec *rec1, *rec2;
 631
 632         if (!up1 || !up2)
 633                 goto error;
 634
 635         if (isl_upoly_is_nan(up1)) {
 636                 isl_upoly_free(up2);
 637                 return up1;
 638         }
 639
 640         if (isl_upoly_is_nan(up2)) {
 641                 isl_upoly_free(up1);
 642                 return up2;
 643         }
 644
 645         if (isl_upoly_is_zero(up1)) {
 646                 isl_upoly_free(up1);
 647                 return up2;
 648         }
 649
 650         if (isl_upoly_is_zero(up2)) {
 651                 isl_upoly_free(up2);
 652                 return up1;
 653         }
 654
 655         if (up1->var < up2->var)
 656                 return isl_upoly_sum(up2, up1);
 657
 658         if (up2->var < up1->var) {
 659                 struct isl_upoly_rec *rec;
 660                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 661                         isl_upoly_free(up1);
 662                         return up2;
 663                 }
 664                 up1 = isl_upoly_cow(up1);
 665                 rec = isl_upoly_as_rec(up1);
 666                 if (!rec)
 667                         goto error;
 668                 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
 669                 if (rec->n == 1)
 670                         up1 = replace_by_constant_term(up1);
 671                 return up1;
 672         }
 673
 674         if (isl_upoly_is_cst(up1))
 675                 return isl_upoly_sum_cst(up1, up2);
 676
 677         rec1 = isl_upoly_as_rec(up1);
 678         rec2 = isl_upoly_as_rec(up2);
 679         if (!rec1 || !rec2)
 680                 goto error;
 681
 682         if (rec1->n < rec2->n)
 683                 return isl_upoly_sum(up2, up1);
 684
 685         up1 = isl_upoly_cow(up1);
 686         rec1 = isl_upoly_as_rec(up1);
 687         if (!rec1)
 688                 goto error;
 689
 690         for (i = rec2->n - 1; i >= 0; --i) {
 691                 rec1->p[i] = isl_upoly_sum(rec1->p[i],
 692                                             isl_upoly_copy(rec2->p[i]));
 693                 if (!rec1->p[i])
 694                         goto error;
 695                 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
 696                         isl_upoly_free(rec1->p[i]);
 697                         rec1->n--;
 698                 }
 699         }
 700
 701         if (rec1->n == 0)
 702                 up1 = replace_by_zero(up1);
 703         else if (rec1->n == 1)
 704                 up1 = replace_by_constant_term(up1);
 705
 706         isl_upoly_free(up2);
 707
 708         return up1;
 709 error:
 710         isl_upoly_free(up1);
 711         isl_upoly_free(up2);
 712         return NULL;
 713 }
 714
 715 __isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
 716         __isl_take struct isl_upoly *up, isl_int v)
 717 {
 718         struct isl_upoly_cst *cst;
 719
 720         up = isl_upoly_cow(up);
 721         if (!up)
 722                 return NULL;
 723
 724         cst = isl_upoly_as_cst(up);
 725
 726         isl_int_addmul(cst->n, cst->d, v);
 727
 728         return up;
 729 }
 730
 731 __isl_give struct isl_upoly *isl_upoly_add_isl_int(
 732         __isl_take struct isl_upoly *up, isl_int v)
 733 {
 734         struct isl_upoly_rec *rec;
 735
 736         if (!up)
 737                 return NULL;
 738
 739         if (isl_upoly_is_cst(up))
 740                 return isl_upoly_cst_add_isl_int(up, v);
 741
 742         up = isl_upoly_cow(up);
 743         rec = isl_upoly_as_rec(up);
 744         if (!rec)
 745                 goto error;
 746
 747         rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
 748         if (!rec->p[0])
 749                 goto error;
 750
 751         return up;
 752 error:
 753         isl_upoly_free(up);
 754         return NULL;
 755 }
 756
 757 __isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
 758         __isl_take struct isl_upoly *up, isl_int v)
 759 {
 760         struct isl_upoly_cst *cst;
 761
 762         if (isl_upoly_is_zero(up))
 763                 return up;
 764
 765         up = isl_upoly_cow(up);
 766         if (!up)
 767                 return NULL;
 768
 769         cst = isl_upoly_as_cst(up);
 770
 771         isl_int_mul(cst->n, cst->n, v);
 772
 773         return up;
 774 }
 775
 776 __isl_give struct isl_upoly *isl_upoly_mul_isl_int(
 777         __isl_take struct isl_upoly *up, isl_int v)
 778 {
 779         int i;
 780         struct isl_upoly_rec *rec;
 781
 782         if (!up)
 783                 return NULL;
 784
 785         if (isl_upoly_is_cst(up))
 786                 return isl_upoly_cst_mul_isl_int(up, v);
 787
 788         up = isl_upoly_cow(up);
 789         rec = isl_upoly_as_rec(up);
 790         if (!rec)
 791                 goto error;
 792
 793         for (i = 0; i < rec->n; ++i) {
 794                 rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
 795                 if (!rec->p[i])
 796                         goto error;
 797         }
 798
 799         return up;
 800 error:
 801         isl_upoly_free(up);
 802         return NULL;
 803 }
 804
 805 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
 806         __isl_take struct isl_upoly *up2)
 807 {
 808         struct isl_upoly_cst *cst1;
 809         struct isl_upoly_cst *cst2;
 810
 811         up1 = isl_upoly_cow(up1);
 812         if (!up1 || !up2)
 813                 goto error;
 814
 815         cst1 = isl_upoly_as_cst(up1);
 816         cst2 = isl_upoly_as_cst(up2);
 817
 818         isl_int_mul(cst1->n, cst1->n, cst2->n);
 819         isl_int_mul(cst1->d, cst1->d, cst2->d);
 820
 821         isl_upoly_cst_reduce(cst1);
 822
 823         isl_upoly_free(up2);
 824         return up1;
 825 error:
 826         isl_upoly_free(up1);
 827         isl_upoly_free(up2);
 828         return NULL;
 829 }
 830
 831 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
 832         __isl_take struct isl_upoly *up2)
 833 {
 834         struct isl_upoly_rec *rec1;
 835         struct isl_upoly_rec *rec2;
 836         struct isl_upoly_rec *res = NULL;
 837         int i, j;
 838         int size;
 839
 840         rec1 = isl_upoly_as_rec(up1);
 841         rec2 = isl_upoly_as_rec(up2);
 842         if (!rec1 || !rec2)
 843                 goto error;
 844         size = rec1->n + rec2->n - 1;
 845         res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
 846         if (!res)
 847                 goto error;
 848
 849         for (i = 0; i < rec1->n; ++i) {
 850                 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
 851                                             isl_upoly_copy(rec1->p[i]));
 852                 if (!res->p[i])
 853                         goto error;
 854                 res->n++;
 855         }
 856         for (; i < size; ++i) {
 857                 res->p[i] = isl_upoly_zero(up1->ctx);
 858                 if (!res->p[i])
 859                         goto error;
 860                 res->n++;
 861         }
 862         for (i = 0; i < rec1->n; ++i) {
 863                 for (j = 1; j < rec2->n; ++j) {
 864                         struct isl_upoly *up;
 865                         up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
 866                                             isl_upoly_copy(rec1->p[i]));
 867                         res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
 868                         if (!res->p[i + j])
 869                                 goto error;
 870                 }
 871         }
 872
 873         isl_upoly_free(up1);
 874         isl_upoly_free(up2);
 875
 876         return &res->up;
 877 error:
 878         isl_upoly_free(up1);
 879         isl_upoly_free(up2);
 880         isl_upoly_free(&res->up);
 881         return NULL;
 882 }
 883
 884 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
 885         __isl_take struct isl_upoly *up2)
 886 {
 887         if (!up1 || !up2)
 888                 goto error;
 889
 890         if (isl_upoly_is_nan(up1)) {
 891                 isl_upoly_free(up2);
 892                 return up1;
 893         }
 894
 895         if (isl_upoly_is_nan(up2)) {
 896                 isl_upoly_free(up1);
 897                 return up2;
 898         }
 899
 900         if (isl_upoly_is_zero(up1)) {
 901                 isl_upoly_free(up2);
 902                 return up1;
 903         }
 904
 905         if (isl_upoly_is_zero(up2)) {
 906                 isl_upoly_free(up1);
 907                 return up2;
 908         }
 909
 910         if (isl_upoly_is_one(up1)) {
 911                 isl_upoly_free(up1);
 912                 return up2;
 913         }
 914
 915         if (isl_upoly_is_one(up2)) {
 916                 isl_upoly_free(up2);
 917                 return up1;
 918         }
 919
 920         if (up1->var < up2->var)
 921                 return isl_upoly_mul(up2, up1);
 922
 923         if (up2->var < up1->var) {
 924                 int i;
 925                 struct isl_upoly_rec *rec;
 926                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 927                         isl_ctx *ctx = up1->ctx;
 928                         isl_upoly_free(up1);
 929                         isl_upoly_free(up2);
 930                         return isl_upoly_nan(ctx);
 931                 }
 932                 up1 = isl_upoly_cow(up1);
 933                 rec = isl_upoly_as_rec(up1);
 934                 if (!rec)
 935                         goto error;
 936
 937                 for (i = 0; i < rec->n; ++i) {
 938                         rec->p[i] = isl_upoly_mul(rec->p[i],
 939                                                     isl_upoly_copy(up2));
 940                         if (!rec->p[i])
 941                                 goto error;
 942                 }
 943                 isl_upoly_free(up2);
 944                 return up1;
 945         }
 946
 947         if (isl_upoly_is_cst(up1))
 948                 return isl_upoly_mul_cst(up1, up2);
 949
 950         return isl_upoly_mul_rec(up1, up2);
 951 error:
 952         isl_upoly_free(up1);
 953         isl_upoly_free(up2);
 954         return NULL;
 955 }
 956
 957 __isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
 958         unsigned power)
 959 {
 960         struct isl_upoly *res;
 961
 962         if (!up)
 963                 return NULL;
 964         if (power == 1)
 965                 return up;
 966
 967         if (power % 2)
 968                 res = isl_upoly_copy(up);
 969         else
 970                 res = isl_upoly_one(up->ctx);
 971
 972         while (power >>= 1) {
 973                 up = isl_upoly_mul(up, isl_upoly_copy(up));
 974                 if (power % 2)
 975                         res = isl_upoly_mul(res, isl_upoly_copy(up));
 976         }
 977
 978         isl_upoly_free(up);
 979         return res;
 980 }
 981
 982 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
 983         unsigned n_div, __isl_take struct isl_upoly *up)
 984 {
 985         struct isl_qpolynomial *qp = NULL;
 986         unsigned total;
 987
 988         if (!dim || !up)
 989                 goto error;
 990
 991         if (!isl_space_is_set(dim))
 992                 isl_die(isl_space_get_ctx(dim), isl_error_invalid,
 993                         "domain of polynomial should be a set", goto error);
 994
 995         total = isl_space_dim(dim, isl_dim_all);
 996
 997         qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
 998         if (!qp)
 999                 goto error;
1000
1001         qp->ref = 1;
1002         qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
1003         if (!qp->div)
1004                 goto error;
1005
1006         qp->dim = dim;
1007         qp->upoly = up;
1008
1009         return qp;
1010 error:
1011         isl_space_free(dim);
1012         isl_upoly_free(up);
1013         isl_qpolynomial_free(qp);
1014         return NULL;
1015 }
1016
1017 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1018 {
1019         if (!qp)
1020                 return NULL;
1021
1022         qp->ref++;
1023         return qp;
1024 }
1025
1026 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1027 {
1028         struct isl_qpolynomial *dup;
1029
1030         if (!qp)
1031                 return NULL;
1032
1033         dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1034                                     isl_upoly_copy(qp->upoly));
1035         if (!dup)
1036                 return NULL;
1037         isl_mat_free(dup->div);
1038         dup->div = isl_mat_copy(qp->div);
1039         if (!dup->div)
1040                 goto error;
1041
1042         return dup;
1043 error:
1044         isl_qpolynomial_free(dup);
1045         return NULL;
1046 }
1047
1048 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1049 {
1050         if (!qp)
1051                 return NULL;
1052
1053         if (qp->ref == 1)
1054                 return qp;
1055         qp->ref--;
1056         return isl_qpolynomial_dup(qp);
1057 }
1058
1059 void *isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
1060 {
1061         if (!qp)
1062                 return NULL;
1063
1064         if (--qp->ref > 0)
1065                 return NULL;
1066
1067         isl_space_free(qp->dim);
1068         isl_mat_free(qp->div);
1069         isl_upoly_free(qp->upoly);
1070
1071         free(qp);
1072         return NULL;
1073 }
1074
1075 __isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
1076 {
1077         int i;
1078         struct isl_upoly_rec *rec;
1079         struct isl_upoly_cst *cst;
1080
1081         rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1082         if (!rec)
1083                 return NULL;
1084         for (i = 0; i < 1 + power; ++i) {
1085                 rec->p[i] = isl_upoly_zero(ctx);
1086                 if (!rec->p[i])
1087                         goto error;
1088                 rec->n++;
1089         }
1090         cst = isl_upoly_as_cst(rec->p[power]);
1091         isl_int_set_si(cst->n, 1);
1092
1093         return &rec->up;
1094 error:
1095         isl_upoly_free(&rec->up);
1096         return NULL;
1097 }
1098
1099 /* r array maps original positions to new positions.
1100  */
1101 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
1102         int *r)
1103 {
1104         int i;
1105         struct isl_upoly_rec *rec;
1106         struct isl_upoly *base;
1107         struct isl_upoly *res;
1108
1109         if (isl_upoly_is_cst(up))
1110                 return up;
1111
1112         rec = isl_upoly_as_rec(up);
1113         if (!rec)
1114                 goto error;
1115
1116         isl_assert(up->ctx, rec->n >= 1, goto error);
1117
1118         base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
1119         res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
1120
1121         for (i = rec->n - 2; i >= 0; --i) {
1122                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1123                 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
1124         }
1125
1126         isl_upoly_free(base);
1127         isl_upoly_free(up);
1128
1129         return res;
1130 error:
1131         isl_upoly_free(up);
1132         return NULL;
1133 }
1134
1135 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
1136 {
1137         int n_row, n_col;
1138         int equal;
1139
1140         isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1141                                 div1->n_col >= div2->n_col, return -1);
1142
1143         if (div1->n_row == div2->n_row)
1144                 return isl_mat_is_equal(div1, div2);
1145
1146         n_row = div1->n_row;
1147         n_col = div1->n_col;
1148         div1->n_row = div2->n_row;
1149         div1->n_col = div2->n_col;
1150
1151         equal = isl_mat_is_equal(div1, div2);
1152
1153         div1->n_row = n_row;
1154         div1->n_col = n_col;
1155
1156         return equal;
1157 }
1158
1159 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1160 {
1161         int li, lj;
1162
1163         li = isl_seq_last_non_zero(div->row[i], div->n_col);
1164         lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1165
1166         if (li != lj)
1167                 return li - lj;
1168
1169         return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1170 }
1171
1172 struct isl_div_sort_info {
1173         isl_mat *div;
1174         int      row;
1175 };
1176
1177 static int div_sort_cmp(const void *p1, const void *p2)
1178 {
1179         const struct isl_div_sort_info *i1, *i2;
1180         i1 = (const struct isl_div_sort_info *) p1;
1181         i2 = (const struct isl_div_sort_info *) p2;
1182
1183         return cmp_row(i1->div, i1->row, i2->row);
1184 }
1185
1186 /* Sort divs and remove duplicates.
1187  */
1188 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1189 {
1190         int i;
1191         int skip;
1192         int len;
1193         struct isl_div_sort_info *array = NULL;
1194         int *pos = NULL, *at = NULL;
1195         int *reordering = NULL;
1196         unsigned div_pos;
1197
1198         if (!qp)
1199                 return NULL;
1200         if (qp->div->n_row <= 1)
1201                 return qp;
1202
1203         div_pos = isl_space_dim(qp->dim, isl_dim_all);
1204
1205         array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1206                                 qp->div->n_row);
1207         pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1208         at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1209         len = qp->div->n_col - 2;
1210         reordering = isl_alloc_array(qp->div->ctx, int, len);
1211         if (!array || !pos || !at || !reordering)
1212                 goto error;
1213
1214         for (i = 0; i < qp->div->n_row; ++i) {
1215                 array[i].div = qp->div;
1216                 array[i].row = i;
1217                 pos[i] = i;
1218                 at[i] = i;
1219         }
1220
1221         qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1222                 div_sort_cmp);
1223
1224         for (i = 0; i < div_pos; ++i)
1225                 reordering[i] = i;
1226
1227         for (i = 0; i < qp->div->n_row; ++i) {
1228                 if (pos[array[i].row] == i)
1229                         continue;
1230                 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1231                 pos[at[i]] = pos[array[i].row];
1232                 at[pos[array[i].row]] = at[i];
1233                 at[i] = array[i].row;
1234                 pos[array[i].row] = i;
1235         }
1236
1237         skip = 0;
1238         for (i = 0; i < len - div_pos; ++i) {
1239                 if (i > 0 &&
1240                     isl_seq_eq(qp->div->row[i - skip - 1],
1241                                qp->div->row[i - skip], qp->div->n_col)) {
1242                         qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1243                         isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1244                                                  2 + div_pos + i - skip);
1245                         qp->div = isl_mat_drop_cols(qp->div,
1246                                                     2 + div_pos + i - skip, 1);
1247                         skip++;
1248                 }
1249                 reordering[div_pos + array[i].row] = div_pos + i - skip;
1250         }
1251
1252         qp->upoly = reorder(qp->upoly, reordering);
1253
1254         if (!qp->upoly || !qp->div)
1255                 goto error;
1256
1257         free(at);
1258         free(pos);
1259         free(array);
1260         free(reordering);
1261
1262         return qp;
1263 error:
1264         free(at);
1265         free(pos);
1266         free(array);
1267         free(reordering);
1268         isl_qpolynomial_free(qp);
1269         return NULL;
1270 }
1271
1272 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1273         int *exp, int first)
1274 {
1275         int i;
1276         struct isl_upoly_rec *rec;
1277
1278         if (isl_upoly_is_cst(up))
1279                 return up;
1280
1281         if (up->var < first)
1282                 return up;
1283
1284         if (exp[up->var - first] == up->var - first)
1285                 return up;
1286
1287         up = isl_upoly_cow(up);
1288         if (!up)
1289                 goto error;
1290
1291         up->var = exp[up->var - first] + first;
1292
1293         rec = isl_upoly_as_rec(up);
1294         if (!rec)
1295                 goto error;
1296
1297         for (i = 0; i < rec->n; ++i) {
1298                 rec->p[i] = expand(rec->p[i], exp, first);
1299                 if (!rec->p[i])
1300                         goto error;
1301         }
1302
1303         return up;
1304 error:
1305         isl_upoly_free(up);
1306         return NULL;
1307 }
1308
1309 static __isl_give isl_qpolynomial *with_merged_divs(
1310         __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1311                                           __isl_take isl_qpolynomial *qp2),
1312         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1313 {
1314         int *exp1 = NULL;
1315         int *exp2 = NULL;
1316         isl_mat *div = NULL;
1317
1318         qp1 = isl_qpolynomial_cow(qp1);
1319         qp2 = isl_qpolynomial_cow(qp2);
1320
1321         if (!qp1 || !qp2)
1322                 goto error;
1323
1324         isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1325                                 qp1->div->n_col >= qp2->div->n_col, goto error);
1326
1327         exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1328         exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1329         if (!exp1 || !exp2)
1330                 goto error;
1331
1332         div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1333         if (!div)
1334                 goto error;
1335
1336         isl_mat_free(qp1->div);
1337         qp1->div = isl_mat_copy(div);
1338         isl_mat_free(qp2->div);
1339         qp2->div = isl_mat_copy(div);
1340
1341         qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1342         qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1343
1344         if (!qp1->upoly || !qp2->upoly)
1345                 goto error;
1346
1347         isl_mat_free(div);
1348         free(exp1);
1349         free(exp2);
1350
1351         return fn(qp1, qp2);
1352 error:
1353         isl_mat_free(div);
1354         free(exp1);
1355         free(exp2);
1356         isl_qpolynomial_free(qp1);
1357         isl_qpolynomial_free(qp2);
1358         return NULL;
1359 }
1360
1361 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1362         __isl_take isl_qpolynomial *qp2)
1363 {
1364         qp1 = isl_qpolynomial_cow(qp1);
1365
1366         if (!qp1 || !qp2)
1367                 goto error;
1368
1369         if (qp1->div->n_row < qp2->div->n_row)
1370                 return isl_qpolynomial_add(qp2, qp1);
1371
1372         isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1373         if (!compatible_divs(qp1->div, qp2->div))
1374                 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1375
1376         qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1377         if (!qp1->upoly)
1378                 goto error;
1379
1380         isl_qpolynomial_free(qp2);
1381
1382         return qp1;
1383 error:
1384         isl_qpolynomial_free(qp1);
1385         isl_qpolynomial_free(qp2);
1386         return NULL;
1387 }
1388
1389 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1390         __isl_keep isl_set *dom,
1391         __isl_take isl_qpolynomial *qp1,
1392         __isl_take isl_qpolynomial *qp2)
1393 {
1394         qp1 = isl_qpolynomial_add(qp1, qp2);
1395         qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1396         return qp1;
1397 }
1398
1399 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1400         __isl_take isl_qpolynomial *qp2)
1401 {
1402         return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1403 }
1404
1405 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1406         __isl_take isl_qpolynomial *qp, isl_int v)
1407 {
1408         if (isl_int_is_zero(v))
1409                 return qp;
1410
1411         qp = isl_qpolynomial_cow(qp);
1412         if (!qp)
1413                 return NULL;
1414
1415         qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1416         if (!qp->upoly)
1417                 goto error;
1418
1419         return qp;
1420 error:
1421         isl_qpolynomial_free(qp);
1422         return NULL;
1423
1424 }
1425
1426 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1427 {
1428         if (!qp)
1429                 return NULL;
1430
1431         return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1432 }
1433
1434 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1435         __isl_take isl_qpolynomial *qp, isl_int v)
1436 {
1437         if (isl_int_is_one(v))
1438                 return qp;
1439
1440         if (qp && isl_int_is_zero(v)) {
1441                 isl_qpolynomial *zero;
1442                 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1443                 isl_qpolynomial_free(qp);
1444                 return zero;
1445         }
1446
1447         qp = isl_qpolynomial_cow(qp);
1448         if (!qp)
1449                 return NULL;
1450
1451         qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1452         if (!qp->upoly)
1453                 goto error;
1454
1455         return qp;
1456 error:
1457         isl_qpolynomial_free(qp);
1458         return NULL;
1459 }
1460
1461 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1462         __isl_take isl_qpolynomial *qp, isl_int v)
1463 {
1464         return isl_qpolynomial_mul_isl_int(qp, v);
1465 }
1466
1467 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1468         __isl_take isl_qpolynomial *qp2)
1469 {
1470         qp1 = isl_qpolynomial_cow(qp1);
1471
1472         if (!qp1 || !qp2)
1473                 goto error;
1474
1475         if (qp1->div->n_row < qp2->div->n_row)
1476                 return isl_qpolynomial_mul(qp2, qp1);
1477
1478         isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1479         if (!compatible_divs(qp1->div, qp2->div))
1480                 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1481
1482         qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1483         if (!qp1->upoly)
1484                 goto error;
1485
1486         isl_qpolynomial_free(qp2);
1487
1488         return qp1;
1489 error:
1490         isl_qpolynomial_free(qp1);
1491         isl_qpolynomial_free(qp2);
1492         return NULL;
1493 }
1494
1495 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1496         unsigned power)
1497 {
1498         qp = isl_qpolynomial_cow(qp);
1499
1500         if (!qp)
1501                 return NULL;
1502
1503         qp->upoly = isl_upoly_pow(qp->upoly, power);
1504         if (!qp->upoly)
1505                 goto error;
1506
1507         return qp;
1508 error:
1509         isl_qpolynomial_free(qp);
1510         return NULL;
1511 }
1512
1513 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1514         __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1515 {
1516         int i;
1517
1518         if (power == 1)
1519                 return pwqp;
1520
1521         pwqp = isl_pw_qpolynomial_cow(pwqp);
1522         if (!pwqp)
1523                 return NULL;
1524
1525         for (i = 0; i < pwqp->n; ++i) {
1526                 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1527                 if (!pwqp->p[i].qp)
1528                         return isl_pw_qpolynomial_free(pwqp);
1529         }
1530
1531         return pwqp;
1532 }
1533
1534 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1535         __isl_take isl_space *dim)
1536 {
1537         if (!dim)
1538                 return NULL;
1539         return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1540 }
1541
1542 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1543         __isl_take isl_space *dim)
1544 {
1545         if (!dim)
1546                 return NULL;
1547         return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1548 }
1549
1550 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1551         __isl_take isl_space *dim)
1552 {
1553         if (!dim)
1554                 return NULL;
1555         return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1556 }
1557
1558 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1559         __isl_take isl_space *dim)
1560 {
1561         if (!dim)
1562                 return NULL;
1563         return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1564 }
1565
1566 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1567         __isl_take isl_space *dim)
1568 {
1569         if (!dim)
1570                 return NULL;
1571         return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1572 }
1573
1574 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1575         __isl_take isl_space *dim,
1576         isl_int v)
1577 {
1578         struct isl_qpolynomial *qp;
1579         struct isl_upoly_cst *cst;
1580
1581         if (!dim)
1582                 return NULL;
1583
1584         qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1585         if (!qp)
1586                 return NULL;
1587
1588         cst = isl_upoly_as_cst(qp->upoly);
1589         isl_int_set(cst->n, v);
1590
1591         return qp;
1592 }
1593
1594 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1595         isl_int *n, isl_int *d)
1596 {
1597         struct isl_upoly_cst *cst;
1598
1599         if (!qp)
1600                 return -1;
1601
1602         if (!isl_upoly_is_cst(qp->upoly))
1603                 return 0;
1604
1605         cst = isl_upoly_as_cst(qp->upoly);
1606         if (!cst)
1607                 return -1;
1608
1609         if (n)
1610                 isl_int_set(*n, cst->n);
1611         if (d)
1612                 isl_int_set(*d, cst->d);
1613
1614         return 1;
1615 }
1616
1617 /* Return the constant term of "up".
1618  */
1619 static __isl_give isl_val *isl_upoly_get_constant_val(
1620         __isl_keep struct isl_upoly *up)
1621 {
1622         struct isl_upoly_cst *cst;
1623
1624         if (!up)
1625                 return NULL;
1626
1627         while (!isl_upoly_is_cst(up)) {
1628                 struct isl_upoly_rec *rec;
1629
1630                 rec = isl_upoly_as_rec(up);
1631                 if (!rec)
1632                         return NULL;
1633                 up = rec->p[0];
1634         }
1635
1636         cst = isl_upoly_as_cst(up);
1637         if (!cst)
1638                 return NULL;
1639         return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1640 }
1641
1642 /* Return the constant term of "qp".
1643  */
1644 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1645         __isl_keep isl_qpolynomial *qp)
1646 {
1647         if (!qp)
1648                 return NULL;
1649
1650         return isl_upoly_get_constant_val(qp->upoly);
1651 }
1652
1653 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1654 {
1655         int is_cst;
1656         struct isl_upoly_rec *rec;
1657
1658         if (!up)
1659                 return -1;
1660
1661         if (up->var < 0)
1662                 return 1;
1663
1664         rec = isl_upoly_as_rec(up);
1665         if (!rec)
1666                 return -1;
1667
1668         if (rec->n > 2)
1669                 return 0;
1670
1671         isl_assert(up->ctx, rec->n > 1, return -1);
1672
1673         is_cst = isl_upoly_is_cst(rec->p[1]);
1674         if (is_cst < 0)
1675                 return -1;
1676         if (!is_cst)
1677                 return 0;
1678
1679         return isl_upoly_is_affine(rec->p[0]);
1680 }
1681
1682 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1683 {
1684         if (!qp)
1685                 return -1;
1686
1687         if (qp->div->n_row > 0)
1688                 return 0;
1689
1690         return isl_upoly_is_affine(qp->upoly);
1691 }
1692
1693 static void update_coeff(__isl_keep isl_vec *aff,
1694         __isl_keep struct isl_upoly_cst *cst, int pos)
1695 {
1696         isl_int gcd;
1697         isl_int f;
1698
1699         if (isl_int_is_zero(cst->n))
1700                 return;
1701
1702         isl_int_init(gcd);
1703         isl_int_init(f);
1704         isl_int_gcd(gcd, cst->d, aff->el[0]);
1705         isl_int_divexact(f, cst->d, gcd);
1706         isl_int_divexact(gcd, aff->el[0], gcd);
1707         isl_seq_scale(aff->el, aff->el, f, aff->size);
1708         isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1709         isl_int_clear(gcd);
1710         isl_int_clear(f);
1711 }
1712
1713 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1714         __isl_keep isl_vec *aff)
1715 {
1716         struct isl_upoly_cst *cst;
1717         struct isl_upoly_rec *rec;
1718
1719         if (!up || !aff)
1720                 return -1;
1721
1722         if (up->var < 0) {
1723                 struct isl_upoly_cst *cst;
1724
1725                 cst = isl_upoly_as_cst(up);
1726                 if (!cst)
1727                         return -1;
1728                 update_coeff(aff, cst, 0);
1729                 return 0;
1730         }
1731
1732         rec = isl_upoly_as_rec(up);
1733         if (!rec)
1734                 return -1;
1735         isl_assert(up->ctx, rec->n == 2, return -1);
1736
1737         cst = isl_upoly_as_cst(rec->p[1]);
1738         if (!cst)
1739                 return -1;
1740         update_coeff(aff, cst, 1 + up->var);
1741
1742         return isl_upoly_update_affine(rec->p[0], aff);
1743 }
1744
1745 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1746         __isl_keep isl_qpolynomial *qp)
1747 {
1748         isl_vec *aff;
1749         unsigned d;
1750
1751         if (!qp)
1752                 return NULL;
1753
1754         d = isl_space_dim(qp->dim, isl_dim_all);
1755         aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1756         if (!aff)
1757                 return NULL;
1758
1759         isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1760         isl_int_set_si(aff->el[0], 1);
1761
1762         if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1763                 goto error;
1764
1765         return aff;
1766 error:
1767         isl_vec_free(aff);
1768         return NULL;
1769 }
1770
1771 int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
1772         __isl_keep isl_qpolynomial *qp2)
1773 {
1774         int equal;
1775
1776         if (!qp1 || !qp2)
1777                 return -1;
1778
1779         equal = isl_space_is_equal(qp1->dim, qp2->dim);
1780         if (equal < 0 || !equal)
1781                 return equal;
1782
1783         equal = isl_mat_is_equal(qp1->div, qp2->div);
1784         if (equal < 0 || !equal)
1785                 return equal;
1786
1787         return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1788 }
1789
1790 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1791 {
1792         int i;
1793         struct isl_upoly_rec *rec;
1794
1795         if (isl_upoly_is_cst(up)) {
1796                 struct isl_upoly_cst *cst;
1797                 cst = isl_upoly_as_cst(up);
1798                 if (!cst)
1799                         return;
1800                 isl_int_lcm(*d, *d, cst->d);
1801                 return;
1802         }
1803
1804         rec = isl_upoly_as_rec(up);
1805         if (!rec)
1806                 return;
1807
1808         for (i = 0; i < rec->n; ++i)
1809                 upoly_update_den(rec->p[i], d);
1810 }
1811
1812 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1813 {
1814         isl_int_set_si(*d, 1);
1815         if (!qp)
1816                 return;
1817         upoly_update_den(qp->upoly, d);
1818 }
1819
1820 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
1821         __isl_take isl_space *dim, int pos, int power)
1822 {
1823         struct isl_ctx *ctx;
1824
1825         if (!dim)
1826                 return NULL;
1827
1828         ctx = dim->ctx;
1829
1830         return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
1831 }
1832
1833 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
1834         enum isl_dim_type type, unsigned pos)
1835 {
1836         if (!dim)
1837                 return NULL;
1838
1839         isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
1840         isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
1841
1842         if (type == isl_dim_set)
1843                 pos += isl_space_dim(dim, isl_dim_param);
1844
1845         return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
1846 error:
1847         isl_space_free(dim);
1848         return NULL;
1849 }
1850
1851 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1852         unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1853 {
1854         int i;
1855         struct isl_upoly_rec *rec;
1856         struct isl_upoly *base, *res;
1857
1858         if (!up)
1859                 return NULL;
1860
1861         if (isl_upoly_is_cst(up))
1862                 return up;
1863
1864         if (up->var < first)
1865                 return up;
1866
1867         rec = isl_upoly_as_rec(up);
1868         if (!rec)
1869                 goto error;
1870
1871         isl_assert(up->ctx, rec->n >= 1, goto error);
1872
1873         if (up->var >= first + n)
1874                 base = isl_upoly_var_pow(up->ctx, up->var, 1);
1875         else
1876                 base = isl_upoly_copy(subs[up->var - first]);
1877
1878         res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
1879         for (i = rec->n - 2; i >= 0; --i) {
1880                 struct isl_upoly *t;
1881                 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
1882                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1883                 res = isl_upoly_sum(res, t);
1884         }
1885
1886         isl_upoly_free(base);
1887         isl_upoly_free(up);
1888
1889         return res;
1890 error:
1891         isl_upoly_free(up);
1892         return NULL;
1893 }
1894
1895 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
1896         isl_int denom, unsigned len)
1897 {
1898         int i;
1899         struct isl_upoly *up;
1900
1901         isl_assert(ctx, len >= 1, return NULL);
1902
1903         up = isl_upoly_rat_cst(ctx, f[0], denom);
1904         for (i = 0; i < len - 1; ++i) {
1905                 struct isl_upoly *t;
1906                 struct isl_upoly *c;
1907
1908                 if (isl_int_is_zero(f[1 + i]))
1909                         continue;
1910
1911                 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
1912                 t = isl_upoly_var_pow(ctx, i, 1);
1913                 t = isl_upoly_mul(c, t);
1914                 up = isl_upoly_sum(up, t);
1915         }
1916
1917         return up;
1918 }
1919
1920 /* Remove common factor of non-constant terms and denominator.
1921  */
1922 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
1923 {
1924         isl_ctx *ctx = qp->div->ctx;
1925         unsigned total = qp->div->n_col - 2;
1926
1927         isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
1928         isl_int_gcd(ctx->normalize_gcd,
1929                     ctx->normalize_gcd, qp->div->row[div][0]);
1930         if (isl_int_is_one(ctx->normalize_gcd))
1931                 return;
1932
1933         isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
1934                             ctx->normalize_gcd, total);
1935         isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
1936                             ctx->normalize_gcd);
1937         isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
1938                             ctx->normalize_gcd);
1939 }
1940
1941 /* Replace the integer division identified by "div" by the polynomial "s".
1942  * The integer division is assumed not to appear in the definition
1943  * of any other integer divisions.
1944  */
1945 static __isl_give isl_qpolynomial *substitute_div(
1946         __isl_take isl_qpolynomial *qp,
1947         int div, __isl_take struct isl_upoly *s)
1948 {
1949         int i;
1950         int total;
1951         int *reordering;
1952
1953         if (!qp || !s)
1954                 goto error;
1955
1956         qp = isl_qpolynomial_cow(qp);
1957         if (!qp)
1958                 goto error;
1959
1960         total = isl_space_dim(qp->dim, isl_dim_all);
1961         qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
1962         if (!qp->upoly)
1963                 goto error;
1964
1965         reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
1966         if (!reordering)
1967                 goto error;
1968         for (i = 0; i < total + div; ++i)
1969                 reordering[i] = i;
1970         for (i = total + div + 1; i < total + qp->div->n_row; ++i)
1971                 reordering[i] = i - 1;
1972         qp->div = isl_mat_drop_rows(qp->div, div, 1);
1973         qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
1974         qp->upoly = reorder(qp->upoly, reordering);
1975         free(reordering);
1976
1977         if (!qp->upoly || !qp->div)
1978                 goto error;
1979
1980         isl_upoly_free(s);
1981         return qp;
1982 error:
1983         isl_qpolynomial_free(qp);
1984         isl_upoly_free(s);
1985         return NULL;
1986 }
1987
1988 /* Replace all integer divisions [e/d] that turn out to not actually be integer
1989  * divisions because d is equal to 1 by their definition, i.e., e.
1990  */
1991 static __isl_give isl_qpolynomial *substitute_non_divs(
1992         __isl_take isl_qpolynomial *qp)
1993 {
1994         int i, j;
1995         int total;
1996         struct isl_upoly *s;
1997
1998         if (!qp)
1999                 return NULL;
2000
2001         total = isl_space_dim(qp->dim, isl_dim_all);
2002         for (i = 0; qp && i < qp->div->n_row; ++i) {
2003                 if (!isl_int_is_one(qp->div->row[i][0]))
2004                         continue;
2005                 for (j = i + 1; j < qp->div->n_row; ++j) {
2006                         if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2007                                 continue;
2008                         isl_seq_combine(qp->div->row[j] + 1,
2009                                 qp->div->ctx->one, qp->div->row[j] + 1,
2010                                 qp->div->row[j][2 + total + i],
2011                                 qp->div->row[i] + 1, 1 + total + i);
2012                         isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2013                         normalize_div(qp, j);
2014                 }
2015                 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2016                                         qp->div->row[i][0], qp->div->n_col - 1);
2017                 qp = substitute_div(qp, i, s);
2018                 --i;
2019         }
2020
2021         return qp;
2022 }
2023
2024 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2025  * with d the denominator.  When replacing the coefficient e of x by
2026  * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2027  * inside the division, so we need to add floor(e/d) * x outside.
2028  * That is, we replace q by q' + floor(e/d) * x and we therefore need
2029  * to adjust the coefficient of x in each later div that depends on the
2030  * current div "div" and also in the affine expression "aff"
2031  * (if it too depends on "div").
2032  */
2033 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2034         __isl_keep isl_vec *aff)
2035 {
2036         int i, j;
2037         isl_int v;
2038         unsigned total = qp->div->n_col - qp->div->n_row - 2;
2039
2040         isl_int_init(v);
2041         for (i = 0; i < 1 + total + div; ++i) {
2042                 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2043                     isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2044                         continue;
2045                 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2046                 isl_int_fdiv_r(qp->div->row[div][1 + i],
2047                                 qp->div->row[div][1 + i], qp->div->row[div][0]);
2048                 if (!isl_int_is_zero(aff->el[1 + total + div]))
2049                         isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
2050                 for (j = div + 1; j < qp->div->n_row; ++j) {
2051                         if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2052                                 continue;
2053                         isl_int_addmul(qp->div->row[j][1 + i],
2054                                         v, qp->div->row[j][2 + total + div]);
2055                 }
2056         }
2057         isl_int_clear(v);
2058 }
2059
2060 /* Check if the last non-zero coefficient is bigger that half of the
2061  * denominator.  If so, we will invert the div to further reduce the number
2062  * of distinct divs that may appear.
2063  * If the last non-zero coefficient is exactly half the denominator,
2064  * then we continue looking for earlier coefficients that are bigger
2065  * than half the denominator.
2066  */
2067 static int needs_invert(__isl_keep isl_mat *div, int row)
2068 {
2069         int i;
2070         int cmp;
2071
2072         for (i = div->n_col - 1; i >= 1; --i) {
2073                 if (isl_int_is_zero(div->row[row][i]))
2074                         continue;
2075                 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2076                 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2077                 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2078                 if (cmp)
2079                         return cmp > 0;
2080                 if (i == 1)
2081                         return 1;
2082         }
2083
2084         return 0;
2085 }
2086
2087 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2088  * We only invert the coefficients of e (and the coefficient of q in
2089  * later divs and in "aff").  After calling this function, the
2090  * coefficients of e should be reduced again.
2091  */
2092 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2093         __isl_keep isl_vec *aff)
2094 {
2095         unsigned total = qp->div->n_col - qp->div->n_row - 2;
2096
2097         isl_seq_neg(qp->div->row[div] + 1,
2098                     qp->div->row[div] + 1, qp->div->n_col - 1);
2099         isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2100         isl_int_add(qp->div->row[div][1],
2101                     qp->div->row[div][1], qp->div->row[div][0]);
2102         if (!isl_int_is_zero(aff->el[1 + total + div]))
2103                 isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
2104         isl_mat_col_mul(qp->div, 2 + total + div,
2105                         qp->div->ctx->negone, 2 + total + div);
2106 }
2107
2108 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2109  * in the interval [0, d-1], with d the denominator and such that the
2110  * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2111  *
2112  * After the reduction, some divs may have become redundant or identical,
2113  * so we call substitute_non_divs and sort_divs.  If these functions
2114  * eliminate divs or merge two or more divs into one, the coefficients
2115  * of the enclosing divs may have to be reduced again, so we call
2116  * ourselves recursively if the number of divs decreases.
2117  */
2118 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2119 {
2120         int i;
2121         isl_vec *aff = NULL;
2122         struct isl_upoly *s;
2123         unsigned n_div;
2124
2125         if (!qp)
2126                 return NULL;
2127
2128         aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
2129         aff = isl_vec_clr(aff);
2130         if (!aff)
2131                 goto error;
2132
2133         isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
2134
2135         for (i = 0; i < qp->div->n_row; ++i) {
2136                 normalize_div(qp, i);
2137                 reduce_div(qp, i, aff);
2138                 if (needs_invert(qp->div, i)) {
2139                         invert_div(qp, i, aff);
2140                         reduce_div(qp, i, aff);
2141                 }
2142         }
2143
2144         s = isl_upoly_from_affine(qp->div->ctx, aff->el,
2145                                   qp->div->ctx->one, aff->size);
2146         qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
2147         isl_upoly_free(s);
2148         if (!qp->upoly)
2149                 goto error;
2150
2151         isl_vec_free(aff);
2152
2153         n_div = qp->div->n_row;
2154         qp = substitute_non_divs(qp);
2155         qp = sort_divs(qp);
2156         if (qp && qp->div->n_row < n_div)
2157                 return reduce_divs(qp);
2158
2159         return qp;
2160 error:
2161         isl_qpolynomial_free(qp);
2162         isl_vec_free(aff);
2163         return NULL;
2164 }
2165
2166 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2167         __isl_take isl_space *dim, const isl_int n, const isl_int d)
2168 {
2169         struct isl_qpolynomial *qp;
2170         struct isl_upoly_cst *cst;
2171
2172         if (!dim)
2173                 return NULL;
2174
2175         qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2176         if (!qp)
2177                 return NULL;
2178
2179         cst = isl_upoly_as_cst(qp->upoly);
2180         isl_int_set(cst->n, n);
2181         isl_int_set(cst->d, d);
2182
2183         return qp;
2184 }
2185
2186 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2187 {
2188         struct isl_upoly_rec *rec;
2189         int i;
2190
2191         if (!up)
2192                 return -1;
2193
2194         if (isl_upoly_is_cst(up))
2195                 return 0;
2196
2197         if (up->var < d)
2198                 active[up->var] = 1;
2199
2200         rec = isl_upoly_as_rec(up);
2201         for (i = 0; i < rec->n; ++i)
2202                 if (up_set_active(rec->p[i], active, d) < 0)
2203                         return -1;
2204
2205         return 0;
2206 }
2207
2208 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2209 {
2210         int i, j;
2211         int d = isl_space_dim(qp->dim, isl_dim_all);
2212
2213         if (!qp || !active)
2214                 return -1;
2215
2216         for (i = 0; i < d; ++i)
2217                 for (j = 0; j < qp->div->n_row; ++j) {
2218                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
2219                                 continue;
2220                         active[i] = 1;
2221                         break;
2222                 }
2223
2224         return up_set_active(qp->upoly, active, d);
2225 }
2226
2227 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2228         enum isl_dim_type type, unsigned first, unsigned n)
2229 {
2230         int i;
2231         int *active = NULL;
2232         int involves = 0;
2233
2234         if (!qp)
2235                 return -1;
2236         if (n == 0)
2237                 return 0;
2238
2239         isl_assert(qp->dim->ctx,
2240                     first + n <= isl_qpolynomial_dim(qp, type), return -1);
2241         isl_assert(qp->dim->ctx, type == isl_dim_param ||
2242                                  type == isl_dim_in, return -1);
2243
2244         active = isl_calloc_array(qp->dim->ctx, int,
2245                                         isl_space_dim(qp->dim, isl_dim_all));
2246         if (set_active(qp, active) < 0)
2247                 goto error;
2248
2249         if (type == isl_dim_in)
2250                 first += isl_space_dim(qp->dim, isl_dim_param);
2251         for (i = 0; i < n; ++i)
2252                 if (active[first + i]) {
2253                         involves = 1;
2254                         break;
2255                 }
2256
2257         free(active);
2258
2259         return involves;
2260 error:
2261         free(active);
2262         return -1;
2263 }
2264
2265 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2266  * of the divs that do appear in the quasi-polynomial.
2267  */
2268 static __isl_give isl_qpolynomial *remove_redundant_divs(
2269         __isl_take isl_qpolynomial *qp)
2270 {
2271         int i, j;
2272         int d;
2273         int len;
2274         int skip;
2275         int *active = NULL;
2276         int *reordering = NULL;
2277         int redundant = 0;
2278         int n_div;
2279         isl_ctx *ctx;
2280
2281         if (!qp)
2282                 return NULL;
2283         if (qp->div->n_row == 0)
2284                 return qp;
2285
2286         d = isl_space_dim(qp->dim, isl_dim_all);
2287         len = qp->div->n_col - 2;
2288         ctx = isl_qpolynomial_get_ctx(qp);
2289         active = isl_calloc_array(ctx, int, len);
2290         if (!active)
2291                 goto error;
2292
2293         if (up_set_active(qp->upoly, active, len) < 0)
2294                 goto error;
2295
2296         for (i = qp->div->n_row - 1; i >= 0; --i) {
2297                 if (!active[d + i]) {
2298                         redundant = 1;
2299                         continue;
2300                 }
2301                 for (j = 0; j < i; ++j) {
2302                         if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2303                                 continue;
2304                         active[d + j] = 1;
2305                         break;
2306                 }
2307         }
2308
2309         if (!redundant) {
2310                 free(active);
2311                 return qp;
2312         }
2313
2314         reordering = isl_alloc_array(qp->div->ctx, int, len);
2315         if (!reordering)
2316                 goto error;
2317
2318         for (i = 0; i < d; ++i)
2319                 reordering[i] = i;
2320
2321         skip = 0;
2322         n_div = qp->div->n_row;
2323         for (i = 0; i < n_div; ++i) {
2324                 if (!active[d + i]) {
2325                         qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2326                         qp->div = isl_mat_drop_cols(qp->div,
2327                                                     2 + d + i - skip, 1);
2328                         skip++;
2329                 }
2330                 reordering[d + i] = d + i - skip;
2331         }
2332
2333         qp->upoly = reorder(qp->upoly, reordering);
2334
2335         if (!qp->upoly || !qp->div)
2336                 goto error;
2337
2338         free(active);
2339         free(reordering);
2340
2341         return qp;
2342 error:
2343         free(active);
2344         free(reordering);
2345         isl_qpolynomial_free(qp);
2346         return NULL;
2347 }
2348
2349 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2350         unsigned first, unsigned n)
2351 {
2352         int i;
2353         struct isl_upoly_rec *rec;
2354
2355         if (!up)
2356                 return NULL;
2357         if (n == 0 || up->var < 0 || up->var < first)
2358                 return up;
2359         if (up->var < first + n) {
2360                 up = replace_by_constant_term(up);
2361                 return isl_upoly_drop(up, first, n);
2362         }
2363         up = isl_upoly_cow(up);
2364         if (!up)
2365                 return NULL;
2366         up->var -= n;
2367         rec = isl_upoly_as_rec(up);
2368         if (!rec)
2369                 goto error;
2370
2371         for (i = 0; i < rec->n; ++i) {
2372                 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2373                 if (!rec->p[i])
2374                         goto error;
2375         }
2376
2377         return up;
2378 error:
2379         isl_upoly_free(up);
2380         return NULL;
2381 }
2382
2383 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2384         __isl_take isl_qpolynomial *qp,
2385         enum isl_dim_type type, unsigned pos, const char *s)
2386 {
2387         qp = isl_qpolynomial_cow(qp);
2388         if (!qp)
2389                 return NULL;
2390         qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2391         if (!qp->dim)
2392                 goto error;
2393         return qp;
2394 error:
2395         isl_qpolynomial_free(qp);
2396         return NULL;
2397 }
2398
2399 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2400         __isl_take isl_qpolynomial *qp,
2401         enum isl_dim_type type, unsigned first, unsigned n)
2402 {
2403         if (!qp)
2404                 return NULL;
2405         if (type == isl_dim_out)
2406                 isl_die(qp->dim->ctx, isl_error_invalid,
2407                         "cannot drop output/set dimension",
2408                         goto error);
2409         if (type == isl_dim_in)
2410                 type = isl_dim_set;
2411         if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2412                 return qp;
2413
2414         qp = isl_qpolynomial_cow(qp);
2415         if (!qp)
2416                 return NULL;
2417
2418         isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2419                         goto error);
2420         isl_assert(qp->dim->ctx, type == isl_dim_param ||
2421                                  type == isl_dim_set, goto error);
2422
2423         qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2424         if (!qp->dim)
2425                 goto error;
2426
2427         if (type == isl_dim_set)
2428                 first += isl_space_dim(qp->dim, isl_dim_param);
2429
2430         qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2431         if (!qp->div)
2432                 goto error;
2433
2434         qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2435         if (!qp->upoly)
2436                 goto error;
2437
2438         return qp;
2439 error:
2440         isl_qpolynomial_free(qp);
2441         return NULL;
2442 }
2443
2444 /* Project the domain of the quasi-polynomial onto its parameter space.
2445  * The quasi-polynomial may not involve any of the domain dimensions.
2446  */
2447 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2448         __isl_take isl_qpolynomial *qp)
2449 {
2450         isl_space *space;
2451         unsigned n;
2452         int involves;
2453
2454         n = isl_qpolynomial_dim(qp, isl_dim_in);
2455         involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2456         if (involves < 0)
2457                 return isl_qpolynomial_free(qp);
2458         if (involves)
2459                 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2460                         "polynomial involves some of the domain dimensions",
2461                         return isl_qpolynomial_free(qp));
2462         qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2463         space = isl_qpolynomial_get_domain_space(qp);
2464         space = isl_space_params(space);
2465         qp = isl_qpolynomial_reset_domain_space(qp, space);
2466         return qp;
2467 }
2468
2469 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2470         __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2471 {
2472         int i, j, k;
2473         isl_int denom;
2474         unsigned total;
2475         unsigned n_div;
2476         struct isl_upoly *up;
2477
2478         if (!eq)
2479                 goto error;
2480         if (eq->n_eq == 0) {
2481                 isl_basic_set_free(eq);
2482                 return qp;
2483         }
2484
2485         qp = isl_qpolynomial_cow(qp);
2486         if (!qp)
2487                 goto error;
2488         qp->div = isl_mat_cow(qp->div);
2489         if (!qp->div)
2490                 goto error;
2491
2492         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2493         n_div = eq->n_div;
2494         isl_int_init(denom);
2495         for (i = 0; i < eq->n_eq; ++i) {
2496                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2497                 if (j < 0 || j == 0 || j >= total)
2498                         continue;
2499
2500                 for (k = 0; k < qp->div->n_row; ++k) {
2501                         if (isl_int_is_zero(qp->div->row[k][1 + j]))
2502                                 continue;
2503                         isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2504                                         &qp->div->row[k][0]);
2505                         normalize_div(qp, k);
2506                 }
2507
2508                 if (isl_int_is_pos(eq->eq[i][j]))
2509                         isl_seq_neg(eq->eq[i], eq->eq[i], total);
2510                 isl_int_abs(denom, eq->eq[i][j]);
2511                 isl_int_set_si(eq->eq[i][j], 0);
2512
2513                 up = isl_upoly_from_affine(qp->dim->ctx,
2514                                                    eq->eq[i], denom, total);
2515                 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2516                 isl_upoly_free(up);
2517         }
2518         isl_int_clear(denom);
2519
2520         if (!qp->upoly)
2521                 goto error;
2522
2523         isl_basic_set_free(eq);
2524
2525         qp = substitute_non_divs(qp);
2526         qp = sort_divs(qp);
2527
2528         return qp;
2529 error:
2530         isl_basic_set_free(eq);
2531         isl_qpolynomial_free(qp);
2532         return NULL;
2533 }
2534
2535 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2536  */
2537 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2538         __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2539 {
2540         if (!qp || !eq)
2541                 goto error;
2542         if (qp->div->n_row > 0)
2543                 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2544         return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2545 error:
2546         isl_basic_set_free(eq);
2547         isl_qpolynomial_free(qp);
2548         return NULL;
2549 }
2550
2551 static __isl_give isl_basic_set *add_div_constraints(
2552         __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2553 {
2554         int i;
2555         unsigned total;
2556
2557         if (!bset || !div)
2558                 goto error;
2559
2560         bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2561         if (!bset)
2562                 goto error;
2563         total = isl_basic_set_total_dim(bset);
2564         for (i = 0; i < div->n_row; ++i)
2565                 if (isl_basic_set_add_div_constraints_var(bset,
2566                                     total - div->n_row + i, div->row[i]) < 0)
2567                         goto error;
2568
2569         isl_mat_free(div);
2570         return bset;
2571 error:
2572         isl_mat_free(div);
2573         isl_basic_set_free(bset);
2574         return NULL;
2575 }
2576
2577 /* Look for equalities among the variables shared by context and qp
2578  * and the integer divisions of qp, if any.
2579  * The equalities are then used to eliminate variables and/or integer
2580  * divisions from qp.
2581  */
2582 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2583         __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2584 {
2585         isl_basic_set *aff;
2586
2587         if (!qp)
2588                 goto error;
2589         if (qp->div->n_row > 0) {
2590                 isl_basic_set *bset;
2591                 context = isl_set_add_dims(context, isl_dim_set,
2592                                             qp->div->n_row);
2593                 bset = isl_basic_set_universe(isl_set_get_space(context));
2594                 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2595                 context = isl_set_intersect(context,
2596                                             isl_set_from_basic_set(bset));
2597         }
2598
2599         aff = isl_set_affine_hull(context);
2600         return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2601 error:
2602         isl_qpolynomial_free(qp);
2603         isl_set_free(context);
2604         return NULL;
2605 }
2606
2607 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2608         __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2609 {
2610         isl_space *space = isl_qpolynomial_get_domain_space(qp);
2611         isl_set *dom_context = isl_set_universe(space);
2612         dom_context = isl_set_intersect_params(dom_context, context);
2613         return isl_qpolynomial_gist(qp, dom_context);
2614 }
2615
2616 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2617         __isl_take isl_qpolynomial *qp)
2618 {
2619         isl_set *dom;
2620
2621         if (!qp)
2622                 return NULL;
2623         if (isl_qpolynomial_is_zero(qp)) {
2624                 isl_space *dim = isl_qpolynomial_get_space(qp);
2625                 isl_qpolynomial_free(qp);
2626                 return isl_pw_qpolynomial_zero(dim);
2627         }
2628
2629         dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2630         return isl_pw_qpolynomial_alloc(dom, qp);
2631 }
2632
2633 #undef PW
2634 #define PW isl_pw_qpolynomial
2635 #undef EL
2636 #define EL isl_qpolynomial
2637 #undef EL_IS_ZERO
2638 #define EL_IS_ZERO is_zero
2639 #undef ZERO
2640 #define ZERO zero
2641 #undef IS_ZERO
2642 #define IS_ZERO is_zero
2643 #undef FIELD
2644 #define FIELD qp
2645 #undef DEFAULT_IS_ZERO
2646 #define DEFAULT_IS_ZERO 1
2647
2648 #define NO_PULLBACK
2649
2650 #include <isl_pw_templ.c>
2651
2652 #undef UNION
2653 #define UNION isl_union_pw_qpolynomial
2654 #undef PART
2655 #define PART isl_pw_qpolynomial
2656 #undef PARTS
2657 #define PARTS pw_qpolynomial
2658 #define ALIGN_DOMAIN
2659
2660 #include <isl_union_templ.c>
2661
2662 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2663 {
2664         if (!pwqp)
2665                 return -1;
2666
2667         if (pwqp->n != -1)
2668                 return 0;
2669
2670         if (!isl_set_plain_is_universe(pwqp->p[0].set))
2671                 return 0;
2672
2673         return isl_qpolynomial_is_one(pwqp->p[0].qp);
2674 }
2675
2676 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2677         __isl_take isl_pw_qpolynomial *pwqp1,
2678         __isl_take isl_pw_qpolynomial *pwqp2)
2679 {
2680         return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2681 }
2682
2683 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2684         __isl_take isl_pw_qpolynomial *pwqp1,
2685         __isl_take isl_pw_qpolynomial *pwqp2)
2686 {
2687         int i, j, n;
2688         struct isl_pw_qpolynomial *res;
2689
2690         if (!pwqp1 || !pwqp2)
2691                 goto error;
2692
2693         isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2694                         goto error);
2695
2696         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2697                 isl_pw_qpolynomial_free(pwqp2);
2698                 return pwqp1;
2699         }
2700
2701         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2702                 isl_pw_qpolynomial_free(pwqp1);
2703                 return pwqp2;
2704         }
2705
2706         if (isl_pw_qpolynomial_is_one(pwqp1)) {
2707                 isl_pw_qpolynomial_free(pwqp1);
2708                 return pwqp2;
2709         }
2710
2711         if (isl_pw_qpolynomial_is_one(pwqp2)) {
2712                 isl_pw_qpolynomial_free(pwqp2);
2713                 return pwqp1;
2714         }
2715
2716         n = pwqp1->n * pwqp2->n;
2717         res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2718
2719         for (i = 0; i < pwqp1->n; ++i) {
2720                 for (j = 0; j < pwqp2->n; ++j) {
2721                         struct isl_set *common;
2722                         struct isl_qpolynomial *prod;
2723                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2724                                                 isl_set_copy(pwqp2->p[j].set));
2725                         if (isl_set_plain_is_empty(common)) {
2726                                 isl_set_free(common);
2727                                 continue;
2728                         }
2729
2730                         prod = isl_qpolynomial_mul(
2731                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
2732                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
2733
2734                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
2735                 }
2736         }
2737
2738         isl_pw_qpolynomial_free(pwqp1);
2739         isl_pw_qpolynomial_free(pwqp2);
2740
2741         return res;
2742 error:
2743         isl_pw_qpolynomial_free(pwqp1);
2744         isl_pw_qpolynomial_free(pwqp2);
2745         return NULL;
2746 }
2747
2748 __isl_give struct isl_upoly *isl_upoly_eval(
2749         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
2750 {
2751         int i;
2752         struct isl_upoly_rec *rec;
2753         struct isl_upoly *res;
2754         struct isl_upoly *base;
2755
2756         if (isl_upoly_is_cst(up)) {
2757                 isl_vec_free(vec);
2758                 return up;
2759         }
2760
2761         rec = isl_upoly_as_rec(up);
2762         if (!rec)
2763                 goto error;
2764
2765         isl_assert(up->ctx, rec->n >= 1, goto error);
2766
2767         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
2768
2769         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2770                                 isl_vec_copy(vec));
2771
2772         for (i = rec->n - 2; i >= 0; --i) {
2773                 res = isl_upoly_mul(res, isl_upoly_copy(base));
2774                 res = isl_upoly_sum(res,
2775                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2776                                                             isl_vec_copy(vec)));
2777         }
2778
2779         isl_upoly_free(base);
2780         isl_upoly_free(up);
2781         isl_vec_free(vec);
2782         return res;
2783 error:
2784         isl_upoly_free(up);
2785         isl_vec_free(vec);
2786         return NULL;
2787 }
2788
2789 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
2790         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
2791 {
2792         isl_vec *ext;
2793         struct isl_upoly *up;
2794         isl_space *dim;
2795
2796         if (!qp || !pnt)
2797                 goto error;
2798         isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2799
2800         if (qp->div->n_row == 0)
2801                 ext = isl_vec_copy(pnt->vec);
2802         else {
2803                 int i;
2804                 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2805                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2806                 if (!ext)
2807                         goto error;
2808
2809                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2810                 for (i = 0; i < qp->div->n_row; ++i) {
2811                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2812                                                 1 + dim + i, &ext->el[1+dim+i]);
2813                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2814                                         qp->div->row[i][0]);
2815                 }
2816         }
2817
2818         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2819         if (!up)
2820                 goto error;
2821
2822         dim = isl_space_copy(qp->dim);
2823         isl_qpolynomial_free(qp);
2824         isl_point_free(pnt);
2825
2826         return isl_qpolynomial_alloc(dim, 0, up);
2827 error:
2828         isl_qpolynomial_free(qp);
2829         isl_point_free(pnt);
2830         return NULL;
2831 }
2832
2833 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2834         __isl_keep struct isl_upoly_cst *cst2)
2835 {
2836         int cmp;
2837         isl_int t;
2838         isl_int_init(t);
2839         isl_int_mul(t, cst1->n, cst2->d);
2840         isl_int_submul(t, cst2->n, cst1->d);
2841         cmp = isl_int_sgn(t);
2842         isl_int_clear(t);
2843         return cmp;
2844 }
2845
2846 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial *qp1,
2847         __isl_keep isl_qpolynomial *qp2)
2848 {
2849         struct isl_upoly_cst *cst1, *cst2;
2850
2851         if (!qp1 || !qp2)
2852                 return -1;
2853         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), return -1);
2854         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), return -1);
2855         if (isl_qpolynomial_is_nan(qp1))
2856                 return -1;
2857         if (isl_qpolynomial_is_nan(qp2))
2858                 return -1;
2859         cst1 = isl_upoly_as_cst(qp1->upoly);
2860         cst2 = isl_upoly_as_cst(qp2->upoly);
2861
2862         return isl_upoly_cmp(cst1, cst2) <= 0;
2863 }
2864
2865 __isl_give isl_qpolynomial *isl_qpolynomial_min_cst(
2866         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2867 {
2868         struct isl_upoly_cst *cst1, *cst2;
2869         int cmp;
2870
2871         if (!qp1 || !qp2)
2872                 goto error;
2873         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2874         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2875         cst1 = isl_upoly_as_cst(qp1->upoly);
2876         cst2 = isl_upoly_as_cst(qp2->upoly);
2877         cmp = isl_upoly_cmp(cst1, cst2);
2878
2879         if (cmp <= 0) {
2880                 isl_qpolynomial_free(qp2);
2881         } else {
2882                 isl_qpolynomial_free(qp1);
2883                 qp1 = qp2;
2884         }
2885         return qp1;
2886 error:
2887         isl_qpolynomial_free(qp1);
2888         isl_qpolynomial_free(qp2);
2889         return NULL;
2890 }
2891
2892 __isl_give isl_qpolynomial *isl_qpolynomial_max_cst(
2893         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2894 {
2895         struct isl_upoly_cst *cst1, *cst2;
2896         int cmp;
2897
2898         if (!qp1 || !qp2)
2899                 goto error;
2900         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2901         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2902         cst1 = isl_upoly_as_cst(qp1->upoly);
2903         cst2 = isl_upoly_as_cst(qp2->upoly);
2904         cmp = isl_upoly_cmp(cst1, cst2);
2905
2906         if (cmp >= 0) {
2907                 isl_qpolynomial_free(qp2);
2908         } else {
2909                 isl_qpolynomial_free(qp1);
2910                 qp1 = qp2;
2911         }
2912         return qp1;
2913 error:
2914         isl_qpolynomial_free(qp1);
2915         isl_qpolynomial_free(qp2);
2916         return NULL;
2917 }
2918
2919 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2920         __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2921         unsigned first, unsigned n)
2922 {
2923         unsigned total;
2924         unsigned g_pos;
2925         int *exp;
2926
2927         if (!qp)
2928                 return NULL;
2929         if (type == isl_dim_out)
2930                 isl_die(qp->div->ctx, isl_error_invalid,
2931                         "cannot insert output/set dimensions",
2932                         goto error);
2933         if (type == isl_dim_in)
2934                 type = isl_dim_set;
2935         if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2936                 return qp;
2937
2938         qp = isl_qpolynomial_cow(qp);
2939         if (!qp)
2940                 return NULL;
2941
2942         isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
2943                     goto error);
2944
2945         g_pos = pos(qp->dim, type) + first;
2946
2947         qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
2948         if (!qp->div)
2949                 goto error;
2950
2951         total = qp->div->n_col - 2;
2952         if (total > g_pos) {
2953                 int i;
2954                 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
2955                 if (!exp)
2956                         goto error;
2957                 for (i = 0; i < total - g_pos; ++i)
2958                         exp[i] = i + n;
2959                 qp->upoly = expand(qp->upoly, exp, g_pos);
2960                 free(exp);
2961                 if (!qp->upoly)
2962                         goto error;
2963         }
2964
2965         qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
2966         if (!qp->dim)
2967                 goto error;
2968
2969         return qp;
2970 error:
2971         isl_qpolynomial_free(qp);
2972         return NULL;
2973 }
2974
2975 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
2976         __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
2977 {
2978         unsigned pos;
2979
2980         pos = isl_qpolynomial_dim(qp, type);
2981
2982         return isl_qpolynomial_insert_dims(qp, type, pos, n);
2983 }
2984
2985 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
2986         __isl_take isl_pw_qpolynomial *pwqp,
2987         enum isl_dim_type type, unsigned n)
2988 {
2989         unsigned pos;
2990
2991         pos = isl_pw_qpolynomial_dim(pwqp, type);
2992
2993         return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
2994 }
2995
2996 static int *reordering_move(isl_ctx *ctx,
2997         unsigned len, unsigned dst, unsigned src, unsigned n)
2998 {
2999         int i;
3000         int *reordering;
3001
3002         reordering = isl_alloc_array(ctx, int, len);
3003         if (!reordering)
3004                 return NULL;
3005
3006         if (dst <= src) {
3007                 for (i = 0; i < dst; ++i)
3008                         reordering[i] = i;
3009                 for (i = 0; i < n; ++i)
3010                         reordering[src + i] = dst + i;
3011                 for (i = 0; i < src - dst; ++i)
3012                         reordering[dst + i] = dst + n + i;
3013                 for (i = 0; i < len - src - n; ++i)
3014                         reordering[src + n + i] = src + n + i;
3015         } else {
3016                 for (i = 0; i < src; ++i)
3017                         reordering[i] = i;
3018                 for (i = 0; i < n; ++i)
3019                         reordering[src + i] = dst + i;
3020                 for (i = 0; i < dst - src; ++i)
3021                         reordering[src + n + i] = src + i;
3022                 for (i = 0; i < len - dst - n; ++i)
3023                         reordering[dst + n + i] = dst + n + i;
3024         }
3025
3026         return reordering;
3027 }
3028
3029 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3030         __isl_take isl_qpolynomial *qp,
3031         enum isl_dim_type dst_type, unsigned dst_pos,
3032         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3033 {
3034         unsigned g_dst_pos;
3035         unsigned g_src_pos;
3036         int *reordering;
3037
3038         qp = isl_qpolynomial_cow(qp);
3039         if (!qp)
3040                 return NULL;
3041
3042         if (dst_type == isl_dim_out || src_type == isl_dim_out)
3043                 isl_die(qp->dim->ctx, isl_error_invalid,
3044                         "cannot move output/set dimension",
3045                         goto error);
3046         if (dst_type == isl_dim_in)
3047                 dst_type = isl_dim_set;
3048         if (src_type == isl_dim_in)
3049                 src_type = isl_dim_set;
3050
3051         isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3052                 goto error);
3053
3054         g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3055         g_src_pos = pos(qp->dim, src_type) + src_pos;
3056         if (dst_type > src_type)
3057                 g_dst_pos -= n;
3058
3059         qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3060         if (!qp->div)
3061                 goto error;
3062         qp = sort_divs(qp);
3063         if (!qp)
3064                 goto error;
3065
3066         reordering = reordering_move(qp->dim->ctx,
3067                                 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3068         if (!reordering)
3069                 goto error;
3070
3071         qp->upoly = reorder(qp->upoly, reordering);
3072         free(reordering);
3073         if (!qp->upoly)
3074                 goto error;
3075
3076         qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3077         if (!qp->dim)
3078                 goto error;
3079
3080         return qp;
3081 error:
3082         isl_qpolynomial_free(qp);
3083         return NULL;
3084 }
3085
3086 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3087         isl_int *f, isl_int denom)
3088 {
3089         struct isl_upoly *up;
3090
3091         dim = isl_space_domain(dim);
3092         if (!dim)
3093                 return NULL;
3094
3095         up = isl_upoly_from_affine(dim->ctx, f, denom,
3096                                         1 + isl_space_dim(dim, isl_dim_all));
3097
3098         return isl_qpolynomial_alloc(dim, 0, up);
3099 }
3100
3101 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3102 {
3103         isl_ctx *ctx;
3104         struct isl_upoly *up;
3105         isl_qpolynomial *qp;
3106
3107         if (!aff)
3108                 return NULL;
3109
3110         ctx = isl_aff_get_ctx(aff);
3111         up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3112                                     aff->v->size - 1);
3113
3114         qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3115                                     aff->ls->div->n_row, up);
3116         if (!qp)
3117                 goto error;
3118
3119         isl_mat_free(qp->div);
3120         qp->div = isl_mat_copy(aff->ls->div);
3121         qp->div = isl_mat_cow(qp->div);
3122         if (!qp->div)
3123                 goto error;
3124
3125         isl_aff_free(aff);
3126         qp = reduce_divs(qp);
3127         qp = remove_redundant_divs(qp);
3128         return qp;
3129 error:
3130         isl_aff_free(aff);
3131         return NULL;
3132 }
3133
3134 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3135         __isl_take isl_pw_aff *pwaff)
3136 {
3137         int i;
3138         isl_pw_qpolynomial *pwqp;
3139
3140         if (!pwaff)
3141                 return NULL;
3142
3143         pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3144                                                 pwaff->n);
3145
3146         for (i = 0; i < pwaff->n; ++i) {
3147                 isl_set *dom;
3148                 isl_qpolynomial *qp;
3149
3150                 dom = isl_set_copy(pwaff->p[i].set);
3151                 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3152                 pwqp = isl_pw_qpolynomial_add_piece(pwqp,  dom, qp);
3153         }
3154
3155         isl_pw_aff_free(pwaff);
3156         return pwqp;
3157 }
3158
3159 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3160         __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3161 {
3162         isl_aff *aff;
3163
3164         aff = isl_constraint_get_bound(c, type, pos);
3165         isl_constraint_free(c);
3166         return isl_qpolynomial_from_aff(aff);
3167 }
3168
3169 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3170  * in "qp" by subs[i].
3171  */
3172 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3173         __isl_take isl_qpolynomial *qp,
3174         enum isl_dim_type type, unsigned first, unsigned n,
3175         __isl_keep isl_qpolynomial **subs)
3176 {
3177         int i;
3178         struct isl_upoly **ups;
3179
3180         if (n == 0)
3181                 return qp;
3182
3183         qp = isl_qpolynomial_cow(qp);
3184         if (!qp)
3185                 return NULL;
3186
3187         if (type == isl_dim_out)
3188                 isl_die(qp->dim->ctx, isl_error_invalid,
3189                         "cannot substitute output/set dimension",
3190                         goto error);
3191         if (type == isl_dim_in)
3192                 type = isl_dim_set;
3193
3194         for (i = 0; i < n; ++i)
3195                 if (!subs[i])
3196                         goto error;
3197
3198         isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3199                         goto error);
3200
3201         for (i = 0; i < n; ++i)
3202                 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3203                                 goto error);
3204
3205         isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3206         for (i = 0; i < n; ++i)
3207                 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3208
3209         first += pos(qp->dim, type);
3210
3211         ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3212         if (!ups)
3213                 goto error;
3214         for (i = 0; i < n; ++i)
3215                 ups[i] = subs[i]->upoly;
3216
3217         qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3218
3219         free(ups);
3220
3221         if (!qp->upoly)
3222                 goto error;
3223
3224         return qp;
3225 error:
3226         isl_qpolynomial_free(qp);
3227         return NULL;
3228 }
3229
3230 /* Extend "bset" with extra set dimensions for each integer division
3231  * in "qp" and then call "fn" with the extended bset and the polynomial
3232  * that results from replacing each of the integer divisions by the
3233  * corresponding extra set dimension.
3234  */
3235 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3236         __isl_keep isl_basic_set *bset,
3237         int (*fn)(__isl_take isl_basic_set *bset,
3238                   __isl_take isl_qpolynomial *poly, void *user), void *user)
3239 {
3240         isl_space *dim;
3241         isl_mat *div;
3242         isl_qpolynomial *poly;
3243
3244         if (!qp || !bset)
3245                 goto error;
3246         if (qp->div->n_row == 0)
3247                 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3248                           user);
3249
3250         div = isl_mat_copy(qp->div);
3251         dim = isl_space_copy(qp->dim);
3252         dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3253         poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3254         bset = isl_basic_set_copy(bset);
3255         bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3256         bset = add_div_constraints(bset, div);
3257
3258         return fn(bset, poly, user);
3259 error:
3260         return -1;
3261 }
3262
3263 /* Return total degree in variables first (inclusive) up to last (exclusive).
3264  */
3265 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3266 {
3267         int deg = -1;
3268         int i;
3269         struct isl_upoly_rec *rec;
3270
3271         if (!up)
3272                 return -2;
3273         if (isl_upoly_is_zero(up))
3274                 return -1;
3275         if (isl_upoly_is_cst(up) || up->var < first)
3276                 return 0;
3277
3278         rec = isl_upoly_as_rec(up);
3279         if (!rec)
3280                 return -2;
3281
3282         for (i = 0; i < rec->n; ++i) {
3283                 int d;
3284
3285                 if (isl_upoly_is_zero(rec->p[i]))
3286                         continue;
3287                 d = isl_upoly_degree(rec->p[i], first, last);
3288                 if (up->var < last)
3289                         d += i;
3290                 if (d > deg)
3291                         deg = d;
3292         }
3293
3294         return deg;
3295 }
3296
3297 /* Return total degree in set variables.
3298  */
3299 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3300 {
3301         unsigned ovar;
3302         unsigned nvar;
3303
3304         if (!poly)
3305                 return -2;
3306
3307         ovar = isl_space_offset(poly->dim, isl_dim_set);
3308         nvar = isl_space_dim(poly->dim, isl_dim_set);
3309         return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3310 }
3311
3312 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3313         unsigned pos, int deg)
3314 {
3315         int i;
3316         struct isl_upoly_rec *rec;
3317
3318         if (!up)
3319                 return NULL;
3320
3321         if (isl_upoly_is_cst(up) || up->var < pos) {
3322                 if (deg == 0)
3323                         return isl_upoly_copy(up);
3324                 else
3325                         return isl_upoly_zero(up->ctx);
3326         }
3327
3328         rec = isl_upoly_as_rec(up);
3329         if (!rec)
3330                 return NULL;
3331
3332         if (up->var == pos) {
3333                 if (deg < rec->n)
3334                         return isl_upoly_copy(rec->p[deg]);
3335                 else
3336                         return isl_upoly_zero(up->ctx);
3337         }
3338
3339         up = isl_upoly_copy(up);
3340         up = isl_upoly_cow(up);
3341         rec = isl_upoly_as_rec(up);
3342         if (!rec)
3343                 goto error;
3344
3345         for (i = 0; i < rec->n; ++i) {
3346                 struct isl_upoly *t;
3347                 t = isl_upoly_coeff(rec->p[i], pos, deg);
3348                 if (!t)
3349                         goto error;
3350                 isl_upoly_free(rec->p[i]);
3351                 rec->p[i] = t;
3352         }
3353
3354         return up;
3355 error:
3356         isl_upoly_free(up);
3357         return NULL;
3358 }
3359
3360 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3361  */
3362 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3363         __isl_keep isl_qpolynomial *qp,
3364         enum isl_dim_type type, unsigned t_pos, int deg)
3365 {
3366         unsigned g_pos;
3367         struct isl_upoly *up;
3368         isl_qpolynomial *c;
3369
3370         if (!qp)
3371                 return NULL;
3372
3373         if (type == isl_dim_out)
3374                 isl_die(qp->div->ctx, isl_error_invalid,
3375                         "output/set dimension does not have a coefficient",
3376                         return NULL);
3377         if (type == isl_dim_in)
3378                 type = isl_dim_set;
3379
3380         isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3381                         return NULL);
3382
3383         g_pos = pos(qp->dim, type) + t_pos;
3384         up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3385
3386         c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3387         if (!c)
3388                 return NULL;
3389         isl_mat_free(c->div);
3390         c->div = isl_mat_copy(qp->div);
3391         if (!c->div)
3392                 goto error;
3393         return c;
3394 error:
3395         isl_qpolynomial_free(c);
3396         return NULL;
3397 }
3398
3399 /* Homogenize the polynomial in the variables first (inclusive) up to
3400  * last (exclusive) by inserting powers of variable first.
3401  * Variable first is assumed not to appear in the input.
3402  */
3403 __isl_give struct isl_upoly *isl_upoly_homogenize(
3404         __isl_take struct isl_upoly *up, int deg, int target,
3405         int first, int last)
3406 {
3407         int i;
3408         struct isl_upoly_rec *rec;
3409
3410         if (!up)
3411                 return NULL;
3412         if (isl_upoly_is_zero(up))
3413                 return up;
3414         if (deg == target)
3415                 return up;
3416         if (isl_upoly_is_cst(up) || up->var < first) {
3417                 struct isl_upoly *hom;
3418
3419                 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3420                 if (!hom)
3421                         goto error;
3422                 rec = isl_upoly_as_rec(hom);
3423                 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3424
3425                 return hom;
3426         }
3427
3428         up = isl_upoly_cow(up);
3429         rec = isl_upoly_as_rec(up);
3430         if (!rec)
3431                 goto error;
3432
3433         for (i = 0; i < rec->n; ++i) {
3434                 if (isl_upoly_is_zero(rec->p[i]))
3435                         continue;
3436                 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3437                                 up->var < last ? deg + i : i, target,
3438                                 first, last);
3439                 if (!rec->p[i])
3440                         goto error;
3441         }
3442
3443         return up;
3444 error:
3445         isl_upoly_free(up);
3446         return NULL;
3447 }
3448
3449 /* Homogenize the polynomial in the set variables by introducing
3450  * powers of an extra set variable at position 0.
3451  */
3452 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3453         __isl_take isl_qpolynomial *poly)
3454 {
3455         unsigned ovar;
3456         unsigned nvar;
3457         int deg = isl_qpolynomial_degree(poly);
3458
3459         if (deg < -1)
3460                 goto error;
3461
3462         poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3463         poly = isl_qpolynomial_cow(poly);
3464         if (!poly)
3465                 goto error;
3466
3467         ovar = isl_space_offset(poly->dim, isl_dim_set);
3468         nvar = isl_space_dim(poly->dim, isl_dim_set);
3469         poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3470                                                 ovar, ovar + nvar);
3471         if (!poly->upoly)
3472                 goto error;
3473
3474         return poly;
3475 error:
3476         isl_qpolynomial_free(poly);
3477         return NULL;
3478 }
3479
3480 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3481         __isl_take isl_mat *div)
3482 {
3483         isl_term *term;
3484         int n;
3485
3486         if (!dim || !div)
3487                 goto error;
3488
3489         n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3490
3491         term = isl_calloc(dim->ctx, struct isl_term,
3492                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
3493         if (!term)
3494                 goto error;
3495
3496         term->ref = 1;
3497         term->dim = dim;
3498         term->div = div;
3499         isl_int_init(term->n);
3500         isl_int_init(term->d);
3501
3502         return term;
3503 error:
3504         isl_space_free(dim);
3505         isl_mat_free(div);
3506         return NULL;
3507 }
3508
3509 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3510 {
3511         if (!term)
3512                 return NULL;
3513
3514         term->ref++;
3515         return term;
3516 }
3517
3518 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3519 {
3520         int i;
3521         isl_term *dup;
3522         unsigned total;
3523
3524         if (!term)
3525                 return NULL;
3526
3527         total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3528
3529         dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3530         if (!dup)
3531                 return NULL;
3532
3533         isl_int_set(dup->n, term->n);
3534         isl_int_set(dup->d, term->d);
3535
3536         for (i = 0; i < total; ++i)
3537                 dup->pow[i] = term->pow[i];
3538
3539         return dup;
3540 }
3541
3542 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3543 {
3544         if (!term)
3545                 return NULL;
3546
3547         if (term->ref == 1)
3548                 return term;
3549         term->ref--;
3550         return isl_term_dup(term);
3551 }
3552
3553 void isl_term_free(__isl_take isl_term *term)
3554 {
3555         if (!term)
3556                 return;
3557
3558         if (--term->ref > 0)
3559                 return;
3560
3561         isl_space_free(term->dim);
3562         isl_mat_free(term->div);
3563         isl_int_clear(term->n);
3564         isl_int_clear(term->d);
3565         free(term);
3566 }
3567
3568 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3569 {
3570         if (!term)
3571                 return 0;
3572
3573         switch (type) {
3574         case isl_dim_param:
3575         case isl_dim_in:
3576         case isl_dim_out:       return isl_space_dim(term->dim, type);
3577         case isl_dim_div:       return term->div->n_row;
3578         case isl_dim_all:       return isl_space_dim(term->dim, isl_dim_all) +
3579                                                                 term->div->n_row;
3580         default:                return 0;
3581         }
3582 }
3583
3584 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3585 {
3586         return term ? term->dim->ctx : NULL;
3587 }
3588
3589 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3590 {
3591         if (!term)
3592                 return;
3593         isl_int_set(*n, term->n);
3594 }
3595
3596 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3597 {
3598         if (!term)
3599                 return;
3600         isl_int_set(*d, term->d);
3601 }
3602
3603 /* Return the coefficient of the term "term".
3604  */
3605 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3606 {
3607         if (!term)
3608                 return NULL;
3609
3610         return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3611                                         term->n, term->d);
3612 }
3613
3614 int isl_term_get_exp(__isl_keep isl_term *term,
3615         enum isl_dim_type type, unsigned pos)
3616 {
3617         if (!term)
3618                 return -1;
3619
3620         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3621
3622         if (type >= isl_dim_set)
3623                 pos += isl_space_dim(term->dim, isl_dim_param);
3624         if (type >= isl_dim_div)
3625                 pos += isl_space_dim(term->dim, isl_dim_set);
3626
3627         return term->pow[pos];
3628 }
3629
3630 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3631 {
3632         isl_local_space *ls;
3633         isl_aff *aff;
3634
3635         if (!term)
3636                 return NULL;
3637
3638         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3639                         return NULL);
3640
3641         ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3642                                         isl_mat_copy(term->div));
3643         aff = isl_aff_alloc(ls);
3644         if (!aff)
3645                 return NULL;
3646
3647         isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3648
3649         aff = isl_aff_normalize(aff);
3650
3651         return aff;
3652 }
3653
3654 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3655         int (*fn)(__isl_take isl_term *term, void *user),
3656         __isl_take isl_term *term, void *user)
3657 {
3658         int i;
3659         struct isl_upoly_rec *rec;
3660
3661         if (!up || !term)
3662                 goto error;
3663
3664         if (isl_upoly_is_zero(up))
3665                 return term;
3666
3667         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3668         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3669         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3670
3671         if (isl_upoly_is_cst(up)) {
3672                 struct isl_upoly_cst *cst;
3673                 cst = isl_upoly_as_cst(up);
3674                 if (!cst)
3675                         goto error;
3676                 term = isl_term_cow(term);
3677                 if (!term)
3678                         goto error;
3679                 isl_int_set(term->n, cst->n);
3680                 isl_int_set(term->d, cst->d);
3681                 if (fn(isl_term_copy(term), user) < 0)
3682                         goto error;
3683                 return term;
3684         }
3685
3686         rec = isl_upoly_as_rec(up);
3687         if (!rec)
3688                 goto error;
3689
3690         for (i = 0; i < rec->n; ++i) {
3691                 term = isl_term_cow(term);
3692                 if (!term)
3693                         goto error;
3694                 term->pow[up->var] = i;
3695                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3696                 if (!term)
3697                         goto error;
3698         }
3699         term->pow[up->var] = 0;
3700
3701         return term;
3702 error:
3703         isl_term_free(term);
3704         return NULL;
3705 }
3706
3707 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3708         int (*fn)(__isl_take isl_term *term, void *user), void *user)
3709 {
3710         isl_term *term;
3711
3712         if (!qp)
3713                 return -1;
3714
3715         term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3716         if (!term)
3717                 return -1;
3718
3719         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3720
3721         isl_term_free(term);
3722
3723         return term ? 0 : -1;
3724 }
3725
3726 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3727 {
3728         struct isl_upoly *up;
3729         isl_qpolynomial *qp;
3730         int i, n;
3731
3732         if (!term)
3733                 return NULL;
3734
3735         n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3736
3737         up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3738         for (i = 0; i < n; ++i) {
3739                 if (!term->pow[i])
3740                         continue;
3741                 up = isl_upoly_mul(up,
3742                         isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3743         }
3744
3745         qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3746         if (!qp)
3747                 goto error;
3748         isl_mat_free(qp->div);
3749         qp->div = isl_mat_copy(term->div);
3750         if (!qp->div)
3751                 goto error;
3752
3753         isl_term_free(term);
3754         return qp;
3755 error:
3756         isl_qpolynomial_free(qp);
3757         isl_term_free(term);
3758         return NULL;
3759 }
3760
3761 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3762         __isl_take isl_space *dim)
3763 {
3764         int i;
3765         int extra;
3766         unsigned total;
3767
3768         if (!qp || !dim)
3769                 goto error;
3770
3771         if (isl_space_is_equal(qp->dim, dim)) {
3772                 isl_space_free(dim);
3773                 return qp;
3774         }
3775
3776         qp = isl_qpolynomial_cow(qp);
3777         if (!qp)
3778                 goto error;
3779
3780         extra = isl_space_dim(dim, isl_dim_set) -
3781                         isl_space_dim(qp->dim, isl_dim_set);
3782         total = isl_space_dim(qp->dim, isl_dim_all);
3783         if (qp->div->n_row) {
3784                 int *exp;
3785
3786                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3787                 if (!exp)
3788                         goto error;
3789                 for (i = 0; i < qp->div->n_row; ++i)
3790                         exp[i] = extra + i;
3791                 qp->upoly = expand(qp->upoly, exp, total);
3792                 free(exp);
3793                 if (!qp->upoly)
3794                         goto error;
3795         }
3796         qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3797         if (!qp->div)
3798                 goto error;
3799         for (i = 0; i < qp->div->n_row; ++i)
3800                 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3801
3802         isl_space_free(qp->dim);
3803         qp->dim = dim;
3804
3805         return qp;
3806 error:
3807         isl_space_free(dim);
3808         isl_qpolynomial_free(qp);
3809         return NULL;
3810 }
3811
3812 /* For each parameter or variable that does not appear in qp,
3813  * first eliminate the variable from all constraints and then set it to zero.
3814  */
3815 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3816         __isl_keep isl_qpolynomial *qp)
3817 {
3818         int *active = NULL;
3819         int i;
3820         int d;
3821         unsigned nparam;
3822         unsigned nvar;
3823
3824         if (!set || !qp)
3825                 goto error;
3826
3827         d = isl_space_dim(set->dim, isl_dim_all);
3828         active = isl_calloc_array(set->ctx, int, d);
3829         if (set_active(qp, active) < 0)
3830                 goto error;
3831
3832         for (i = 0; i < d; ++i)
3833                 if (!active[i])
3834                         break;
3835
3836         if (i == d) {
3837                 free(active);
3838                 return set;
3839         }
3840
3841         nparam = isl_space_dim(set->dim, isl_dim_param);
3842         nvar = isl_space_dim(set->dim, isl_dim_set);
3843         for (i = 0; i < nparam; ++i) {
3844                 if (active[i])
3845                         continue;
3846                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3847                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3848         }
3849         for (i = 0; i < nvar; ++i) {
3850                 if (active[nparam + i])
3851                         continue;
3852                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3853                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3854         }
3855
3856         free(active);
3857
3858         return set;
3859 error:
3860         free(active);
3861         isl_set_free(set);
3862         return NULL;
3863 }
3864
3865 struct isl_opt_data {
3866         isl_qpolynomial *qp;
3867         int first;
3868         isl_qpolynomial *opt;
3869         int max;
3870 };
3871
3872 static int opt_fn(__isl_take isl_point *pnt, void *user)
3873 {
3874         struct isl_opt_data *data = (struct isl_opt_data *)user;
3875         isl_qpolynomial *val;
3876
3877         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3878         if (data->first) {
3879                 data->first = 0;
3880                 data->opt = val;
3881         } else if (data->max) {
3882                 data->opt = isl_qpolynomial_max_cst(data->opt, val);
3883         } else {
3884                 data->opt = isl_qpolynomial_min_cst(data->opt, val);
3885         }
3886
3887         return 0;
3888 }
3889
3890 __isl_give isl_qpolynomial *isl_qpolynomial_opt_on_domain(
3891         __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3892 {
3893         struct isl_opt_data data = { NULL, 1, NULL, max };
3894
3895         if (!set || !qp)
3896                 goto error;
3897
3898         if (isl_upoly_is_cst(qp->upoly)) {
3899                 isl_set_free(set);
3900                 return qp;
3901         }
3902
3903         set = fix_inactive(set, qp);
3904
3905         data.qp = qp;
3906         if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3907                 goto error;
3908
3909         if (data.first) {
3910                 isl_space *space = isl_qpolynomial_get_domain_space(qp);
3911                 data.opt = isl_qpolynomial_zero_on_domain(space);
3912         }
3913
3914         isl_set_free(set);
3915         isl_qpolynomial_free(qp);
3916         return data.opt;
3917 error:
3918         isl_set_free(set);
3919         isl_qpolynomial_free(qp);
3920         isl_qpolynomial_free(data.opt);
3921         return NULL;
3922 }
3923
3924 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
3925         __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
3926 {
3927         int i;
3928         int n_sub;
3929         isl_ctx *ctx;
3930         struct isl_upoly **subs;
3931         isl_mat *mat, *diag;
3932
3933         qp = isl_qpolynomial_cow(qp);
3934         if (!qp || !morph)
3935                 goto error;
3936
3937         ctx = qp->dim->ctx;
3938         isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
3939
3940         n_sub = morph->inv->n_row - 1;
3941         if (morph->inv->n_row != morph->inv->n_col)
3942                 n_sub += qp->div->n_row;
3943         subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
3944         if (!subs)
3945                 goto error;
3946
3947         for (i = 0; 1 + i < morph->inv->n_row; ++i)
3948                 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
3949                                         morph->inv->row[0][0], morph->inv->n_col);
3950         if (morph->inv->n_row != morph->inv->n_col)
3951                 for (i = 0; i < qp->div->n_row; ++i)
3952                         subs[morph->inv->n_row - 1 + i] =
3953                             isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
3954
3955         qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
3956
3957         for (i = 0; i < n_sub; ++i)
3958                 isl_upoly_free(subs[i]);
3959         free(subs);
3960
3961         diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
3962         mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
3963         diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
3964         mat = isl_mat_diagonal(mat, diag);
3965         qp->div = isl_mat_product(qp->div, mat);
3966         isl_space_free(qp->dim);
3967         qp->dim = isl_space_copy(morph->ran->dim);
3968
3969         if (!qp->upoly || !qp->div || !qp->dim)
3970                 goto error;
3971
3972         isl_morph_free(morph);
3973
3974         return qp;
3975 error:
3976         isl_qpolynomial_free(qp);
3977         isl_morph_free(morph);
3978         return NULL;
3979 }
3980
3981 static int neg_entry(void **entry, void *user)
3982 {
3983         isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
3984
3985         *pwqp = isl_pw_qpolynomial_neg(*pwqp);
3986
3987         return *pwqp ? 0 : -1;
3988 }
3989
3990 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
3991         __isl_take isl_union_pw_qpolynomial *upwqp)
3992 {
3993         upwqp = isl_union_pw_qpolynomial_cow(upwqp);
3994         if (!upwqp)
3995                 return NULL;
3996
3997         if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
3998                                    &neg_entry, NULL) < 0)
3999                 goto error;
4000
4001         return upwqp;
4002 error:
4003         isl_union_pw_qpolynomial_free(upwqp);
4004         return NULL;
4005 }
4006
4007 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4008         __isl_take isl_union_pw_qpolynomial *upwqp1,
4009         __isl_take isl_union_pw_qpolynomial *upwqp2)
4010 {
4011         return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
4012 }
4013
4014 /* Reorder the columns of the given div definitions according to the
4015  * given reordering.
4016  */
4017 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4018         __isl_take isl_reordering *r)
4019 {
4020         int i, j;
4021         isl_mat *mat;
4022         int extra;
4023
4024         if (!div || !r)
4025                 goto error;
4026
4027         extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4028         mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4029         if (!mat)
4030                 goto error;
4031
4032         for (i = 0; i < div->n_row; ++i) {
4033                 isl_seq_cpy(mat->row[i], div->row[i], 2);
4034                 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4035                 for (j = 0; j < r->len; ++j)
4036                         isl_int_set(mat->row[i][2 + r->pos[j]],
4037                                     div->row[i][2 + j]);
4038         }
4039
4040         isl_reordering_free(r);
4041         isl_mat_free(div);
4042         return mat;
4043 error:
4044         isl_reordering_free(r);
4045         isl_mat_free(div);
4046         return NULL;
4047 }
4048
4049 /* Reorder the dimension of "qp" according to the given reordering.
4050  */
4051 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4052         __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4053 {
4054         qp = isl_qpolynomial_cow(qp);
4055         if (!qp)
4056                 goto error;
4057
4058         r = isl_reordering_extend(r, qp->div->n_row);
4059         if (!r)
4060                 goto error;
4061
4062         qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4063         if (!qp->div)
4064                 goto error;
4065
4066         qp->upoly = reorder(qp->upoly, r->pos);
4067         if (!qp->upoly)
4068                 goto error;
4069
4070         qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4071
4072         isl_reordering_free(r);
4073         return qp;
4074 error:
4075         isl_qpolynomial_free(qp);
4076         isl_reordering_free(r);
4077         return NULL;
4078 }
4079
4080 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4081         __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4082 {
4083         if (!qp || !model)
4084                 goto error;
4085
4086         if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4087                 isl_reordering *exp;
4088
4089                 model = isl_space_drop_dims(model, isl_dim_in,
4090                                         0, isl_space_dim(model, isl_dim_in));
4091                 model = isl_space_drop_dims(model, isl_dim_out,
4092                                         0, isl_space_dim(model, isl_dim_out));
4093                 exp = isl_parameter_alignment_reordering(qp->dim, model);
4094                 exp = isl_reordering_extend_space(exp,
4095                                         isl_qpolynomial_get_domain_space(qp));
4096                 qp = isl_qpolynomial_realign_domain(qp, exp);
4097         }
4098
4099         isl_space_free(model);
4100         return qp;
4101 error:
4102         isl_space_free(model);
4103         isl_qpolynomial_free(qp);
4104         return NULL;
4105 }
4106
4107 struct isl_split_periods_data {
4108         int max_periods;
4109         isl_pw_qpolynomial *res;
4110 };
4111
4112 /* Create a slice where the integer division "div" has the fixed value "v".
4113  * In particular, if "div" refers to floor(f/m), then create a slice
4114  *
4115  *      m v <= f <= m v + (m - 1)
4116  *
4117  * or
4118  *
4119  *      f - m v >= 0
4120  *      -f + m v + (m - 1) >= 0
4121  */
4122 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4123         __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4124 {
4125         int total;
4126         isl_basic_set *bset = NULL;
4127         int k;
4128
4129         if (!dim || !qp)
4130                 goto error;
4131
4132         total = isl_space_dim(dim, isl_dim_all);
4133         bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4134
4135         k = isl_basic_set_alloc_inequality(bset);
4136         if (k < 0)
4137                 goto error;
4138         isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4139         isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4140
4141         k = isl_basic_set_alloc_inequality(bset);
4142         if (k < 0)
4143                 goto error;
4144         isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4145         isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4146         isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4147         isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4148
4149         isl_space_free(dim);
4150         return isl_set_from_basic_set(bset);
4151 error:
4152         isl_basic_set_free(bset);
4153         isl_space_free(dim);
4154         return NULL;
4155 }
4156
4157 static int split_periods(__isl_take isl_set *set,
4158         __isl_take isl_qpolynomial *qp, void *user);
4159
4160 /* Create a slice of the domain "set" such that integer division "div"
4161  * has the fixed value "v" and add the results to data->res,
4162  * replacing the integer division by "v" in "qp".
4163  */
4164 static int set_div(__isl_take isl_set *set,
4165         __isl_take isl_qpolynomial *qp, int div, isl_int v,
4166         struct isl_split_periods_data *data)
4167 {
4168         int i;
4169         int total;
4170         isl_set *slice;
4171         struct isl_upoly *cst;
4172
4173         slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4174         set = isl_set_intersect(set, slice);
4175
4176         if (!qp)
4177                 goto error;
4178
4179         total = isl_space_dim(qp->dim, isl_dim_all);
4180
4181         for (i = div + 1; i < qp->div->n_row; ++i) {
4182                 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4183                         continue;
4184                 isl_int_addmul(qp->div->row[i][1],
4185                                 qp->div->row[i][2 + total + div], v);
4186                 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4187         }
4188
4189         cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4190         qp = substitute_div(qp, div, cst);
4191
4192         return split_periods(set, qp, data);
4193 error:
4194         isl_set_free(set);
4195         isl_qpolynomial_free(qp);
4196         return -1;
4197 }
4198
4199 /* Split the domain "set" such that integer division "div"
4200  * has a fixed value (ranging from "min" to "max") on each slice
4201  * and add the results to data->res.
4202  */
4203 static int split_div(__isl_take isl_set *set,
4204         __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4205         struct isl_split_periods_data *data)
4206 {
4207         for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4208                 isl_set *set_i = isl_set_copy(set);
4209                 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4210
4211                 if (set_div(set_i, qp_i, div, min, data) < 0)
4212                         goto error;
4213         }
4214         isl_set_free(set);
4215         isl_qpolynomial_free(qp);
4216         return 0;
4217 error:
4218         isl_set_free(set);
4219         isl_qpolynomial_free(qp);
4220         return -1;
4221 }
4222
4223 /* If "qp" refers to any integer division
4224  * that can only attain "max_periods" distinct values on "set"
4225  * then split the domain along those distinct values.
4226  * Add the results (or the original if no splitting occurs)
4227  * to data->res.
4228  */
4229 static int split_periods(__isl_take isl_set *set,
4230         __isl_take isl_qpolynomial *qp, void *user)
4231 {
4232         int i;
4233         isl_pw_qpolynomial *pwqp;
4234         struct isl_split_periods_data *data;
4235         isl_int min, max;
4236         int total;
4237         int r = 0;
4238
4239         data = (struct isl_split_periods_data *)user;
4240
4241         if (!set || !qp)
4242                 goto error;
4243
4244         if (qp->div->n_row == 0) {
4245                 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4246                 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4247                 return 0;
4248         }
4249
4250         isl_int_init(min);
4251         isl_int_init(max);
4252         total = isl_space_dim(qp->dim, isl_dim_all);
4253         for (i = 0; i < qp->div->n_row; ++i) {
4254                 enum isl_lp_result lp_res;
4255
4256                 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4257                                                 qp->div->n_row) != -1)
4258                         continue;
4259
4260                 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4261                                           set->ctx->one, &min, NULL, NULL);
4262                 if (lp_res == isl_lp_error)
4263                         goto error2;
4264                 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4265                         continue;
4266                 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4267
4268                 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4269                                           set->ctx->one, &max, NULL, NULL);
4270                 if (lp_res == isl_lp_error)
4271                         goto error2;
4272                 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4273                         continue;
4274                 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4275
4276                 isl_int_sub(max, max, min);
4277                 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4278                         isl_int_add(max, max, min);
4279                         break;
4280                 }
4281         }
4282
4283         if (i < qp->div->n_row) {
4284                 r = split_div(set, qp, i, min, max, data);
4285         } else {
4286                 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4287                 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4288         }
4289
4290         isl_int_clear(max);
4291         isl_int_clear(min);
4292
4293         return r;
4294 error2:
4295         isl_int_clear(max);
4296         isl_int_clear(min);
4297 error:
4298         isl_set_free(set);
4299         isl_qpolynomial_free(qp);
4300         return -1;
4301 }
4302
4303 /* If any quasi-polynomial in pwqp refers to any integer division
4304  * that can only attain "max_periods" distinct values on its domain
4305  * then split the domain along those distinct values.
4306  */
4307 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4308         __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4309 {
4310         struct isl_split_periods_data data;
4311
4312         data.max_periods = max_periods;
4313         data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4314
4315         if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4316                 goto error;
4317
4318         isl_pw_qpolynomial_free(pwqp);
4319
4320         return data.res;
4321 error:
4322         isl_pw_qpolynomial_free(data.res);
4323         isl_pw_qpolynomial_free(pwqp);
4324         return NULL;
4325 }
4326
4327 /* Construct a piecewise quasipolynomial that is constant on the given
4328  * domain.  In particular, it is
4329  *      0       if cst == 0
4330  *      1       if cst == 1
4331  *  infinity    if cst == -1
4332  */
4333 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4334         __isl_take isl_basic_set *bset, int cst)
4335 {
4336         isl_space *dim;
4337         isl_qpolynomial *qp;
4338
4339         if (!bset)
4340                 return NULL;
4341
4342         bset = isl_basic_set_params(bset);
4343         dim = isl_basic_set_get_space(bset);
4344         if (cst < 0)
4345                 qp = isl_qpolynomial_infty_on_domain(dim);
4346         else if (cst == 0)
4347                 qp = isl_qpolynomial_zero_on_domain(dim);
4348         else
4349                 qp = isl_qpolynomial_one_on_domain(dim);
4350         return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4351 }
4352
4353 /* Factor bset, call fn on each of the factors and return the product.
4354  *
4355  * If no factors can be found, simply call fn on the input.
4356  * Otherwise, construct the factors based on the factorizer,
4357  * call fn on each factor and compute the product.
4358  */
4359 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4360         __isl_take isl_basic_set *bset,
4361         __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4362 {
4363         int i, n;
4364         isl_space *dim;
4365         isl_set *set;
4366         isl_factorizer *f;
4367         isl_qpolynomial *qp;
4368         isl_pw_qpolynomial *pwqp;
4369         unsigned nparam;
4370         unsigned nvar;
4371
4372         f = isl_basic_set_factorizer(bset);
4373         if (!f)
4374                 goto error;
4375         if (f->n_group == 0) {
4376                 isl_factorizer_free(f);
4377                 return fn(bset);
4378         }
4379
4380         nparam = isl_basic_set_dim(bset, isl_dim_param);
4381         nvar = isl_basic_set_dim(bset, isl_dim_set);
4382
4383         dim = isl_basic_set_get_space(bset);
4384         dim = isl_space_domain(dim);
4385         set = isl_set_universe(isl_space_copy(dim));
4386         qp = isl_qpolynomial_one_on_domain(dim);
4387         pwqp = isl_pw_qpolynomial_alloc(set, qp);
4388
4389         bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4390
4391         for (i = 0, n = 0; i < f->n_group; ++i) {
4392                 isl_basic_set *bset_i;
4393                 isl_pw_qpolynomial *pwqp_i;
4394
4395                 bset_i = isl_basic_set_copy(bset);
4396                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4397                             nparam + n + f->len[i], nvar - n - f->len[i]);
4398                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4399                             nparam, n);
4400                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4401                             n + f->len[i], nvar - n - f->len[i]);
4402                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4403
4404                 pwqp_i = fn(bset_i);
4405                 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4406
4407                 n += f->len[i];
4408         }
4409
4410         isl_basic_set_free(bset);
4411         isl_factorizer_free(f);
4412
4413         return pwqp;
4414 error:
4415         isl_basic_set_free(bset);
4416         return NULL;
4417 }
4418
4419 /* Factor bset, call fn on each of the factors and return the product.
4420  * The function is assumed to evaluate to zero on empty domains,
4421  * to one on zero-dimensional domains and to infinity on unbounded domains
4422  * and will not be called explicitly on zero-dimensional or unbounded domains.
4423  *
4424  * We first check for some special cases and remove all equalities.
4425  * Then we hand over control to compressed_multiplicative_call.
4426  */
4427 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4428         __isl_take isl_basic_set *bset,
4429         __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4430 {
4431         int bounded;
4432         isl_morph *morph;
4433         isl_pw_qpolynomial *pwqp;
4434
4435         if (!bset)
4436                 return NULL;
4437
4438         if (isl_basic_set_plain_is_empty(bset))
4439                 return constant_on_domain(bset, 0);
4440
4441         if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4442                 return constant_on_domain(bset, 1);
4443
4444         bounded = isl_basic_set_is_bounded(bset);
4445         if (bounded < 0)
4446                 goto error;
4447         if (!bounded)
4448                 return constant_on_domain(bset, -1);
4449
4450         if (bset->n_eq == 0)
4451                 return compressed_multiplicative_call(bset, fn);
4452
4453         morph = isl_basic_set_full_compression(bset);
4454         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4455
4456         pwqp = compressed_multiplicative_call(bset, fn);
4457
4458         morph = isl_morph_dom_params(morph);
4459         morph = isl_morph_ran_params(morph);
4460         morph = isl_morph_inverse(morph);
4461
4462         pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4463
4464         return pwqp;
4465 error:
4466         isl_basic_set_free(bset);
4467         return NULL;
4468 }
4469
4470 /* Drop all floors in "qp", turning each integer division [a/m] into
4471  * a rational division a/m.  If "down" is set, then the integer division
4472  * is replaces by (a-(m-1))/m instead.
4473  */
4474 static __isl_give isl_qpolynomial *qp_drop_floors(
4475         __isl_take isl_qpolynomial *qp, int down)
4476 {
4477         int i;
4478         struct isl_upoly *s;
4479
4480         if (!qp)
4481                 return NULL;
4482         if (qp->div->n_row == 0)
4483                 return qp;
4484
4485         qp = isl_qpolynomial_cow(qp);
4486         if (!qp)
4487                 return NULL;
4488
4489         for (i = qp->div->n_row - 1; i >= 0; --i) {
4490                 if (down) {
4491                         isl_int_sub(qp->div->row[i][1],
4492                                     qp->div->row[i][1], qp->div->row[i][0]);
4493                         isl_int_add_ui(qp->div->row[i][1],
4494                                        qp->div->row[i][1], 1);
4495                 }
4496                 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4497                                         qp->div->row[i][0], qp->div->n_col - 1);
4498                 qp = substitute_div(qp, i, s);
4499                 if (!qp)
4500                         return NULL;
4501         }
4502
4503         return qp;
4504 }
4505
4506 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4507  * a rational division a/m.
4508  */
4509 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4510         __isl_take isl_pw_qpolynomial *pwqp)
4511 {
4512         int i;
4513
4514         if (!pwqp)
4515                 return NULL;
4516
4517         if (isl_pw_qpolynomial_is_zero(pwqp))
4518                 return pwqp;
4519
4520         pwqp = isl_pw_qpolynomial_cow(pwqp);
4521         if (!pwqp)
4522                 return NULL;
4523
4524         for (i = 0; i < pwqp->n; ++i) {
4525                 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4526                 if (!pwqp->p[i].qp)
4527                         goto error;
4528         }
4529
4530         return pwqp;
4531 error:
4532         isl_pw_qpolynomial_free(pwqp);
4533         return NULL;
4534 }
4535
4536 /* Adjust all the integer divisions in "qp" such that they are at least
4537  * one over the given orthant (identified by "signs").  This ensures
4538  * that they will still be non-negative even after subtracting (m-1)/m.
4539  *
4540  * In particular, f is replaced by f' + v, changing f = [a/m]
4541  * to f' = [(a - m v)/m].
4542  * If the constant term k in a is smaller than m,
4543  * the constant term of v is set to floor(k/m) - 1.
4544  * For any other term, if the coefficient c and the variable x have
4545  * the same sign, then no changes are needed.
4546  * Otherwise, if the variable is positive (and c is negative),
4547  * then the coefficient of x in v is set to floor(c/m).
4548  * If the variable is negative (and c is positive),
4549  * then the coefficient of x in v is set to ceil(c/m).
4550  */
4551 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4552         int *signs)
4553 {
4554         int i, j;
4555         int total;
4556         isl_vec *v = NULL;
4557         struct isl_upoly *s;
4558
4559         qp = isl_qpolynomial_cow(qp);
4560         if (!qp)
4561                 return NULL;
4562         qp->div = isl_mat_cow(qp->div);
4563         if (!qp->div)
4564                 goto error;
4565
4566         total = isl_space_dim(qp->dim, isl_dim_all);
4567         v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4568
4569         for (i = 0; i < qp->div->n_row; ++i) {
4570                 isl_int *row = qp->div->row[i];
4571                 v = isl_vec_clr(v);
4572                 if (!v)
4573                         goto error;
4574                 if (isl_int_lt(row[1], row[0])) {
4575                         isl_int_fdiv_q(v->el[0], row[1], row[0]);
4576                         isl_int_sub_ui(v->el[0], v->el[0], 1);
4577                         isl_int_submul(row[1], row[0], v->el[0]);
4578                 }
4579                 for (j = 0; j < total; ++j) {
4580                         if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4581                                 continue;
4582                         if (signs[j] < 0)
4583                                 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4584                         else
4585                                 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4586                         isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4587                 }
4588                 for (j = 0; j < i; ++j) {
4589                         if (isl_int_sgn(row[2 + total + j]) >= 0)
4590                                 continue;
4591                         isl_int_fdiv_q(v->el[1 + total + j],
4592                                         row[2 + total + j], row[0]);
4593                         isl_int_submul(row[2 + total + j],
4594                                         row[0], v->el[1 + total + j]);
4595                 }
4596                 for (j = i + 1; j < qp->div->n_row; ++j) {
4597                         if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4598                                 continue;
4599                         isl_seq_combine(qp->div->row[j] + 1,
4600                                 qp->div->ctx->one, qp->div->row[j] + 1,
4601                                 qp->div->row[j][2 + total + i], v->el, v->size);
4602                 }
4603                 isl_int_set_si(v->el[1 + total + i], 1);
4604                 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4605                                         qp->div->ctx->one, v->size);
4606                 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4607                 isl_upoly_free(s);
4608                 if (!qp->upoly)
4609                         goto error;
4610         }
4611
4612         isl_vec_free(v);
4613         return qp;
4614 error:
4615         isl_vec_free(v);
4616         isl_qpolynomial_free(qp);
4617         return NULL;
4618 }
4619
4620 struct isl_to_poly_data {
4621         int sign;
4622         isl_pw_qpolynomial *res;
4623         isl_qpolynomial *qp;
4624 };
4625
4626 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4627  * We first make all integer divisions positive and then split the
4628  * quasipolynomials into terms with sign data->sign (the direction
4629  * of the requested approximation) and terms with the opposite sign.
4630  * In the first set of terms, each integer division [a/m] is
4631  * overapproximated by a/m, while in the second it is underapproximated
4632  * by (a-(m-1))/m.
4633  */
4634 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4635         void *user)
4636 {
4637         struct isl_to_poly_data *data = user;
4638         isl_pw_qpolynomial *t;
4639         isl_qpolynomial *qp, *up, *down;
4640
4641         qp = isl_qpolynomial_copy(data->qp);
4642         qp = make_divs_pos(qp, signs);
4643
4644         up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4645         up = qp_drop_floors(up, 0);
4646         down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4647         down = qp_drop_floors(down, 1);
4648
4649         isl_qpolynomial_free(qp);
4650         qp = isl_qpolynomial_add(up, down);
4651
4652         t = isl_pw_qpolynomial_alloc(orthant, qp);
4653         data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4654
4655         return 0;
4656 }
4657
4658 /* Approximate each quasipolynomial by a polynomial.  If "sign" is positive,
4659  * the polynomial will be an overapproximation.  If "sign" is negative,
4660  * it will be an underapproximation.  If "sign" is zero, the approximation
4661  * will lie somewhere in between.
4662  *
4663  * In particular, is sign == 0, we simply drop the floors, turning
4664  * the integer divisions into rational divisions.
4665  * Otherwise, we split the domains into orthants, make all integer divisions
4666  * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4667  * depending on the requested sign and the sign of the term in which
4668  * the integer division appears.
4669  */
4670 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4671         __isl_take isl_pw_qpolynomial *pwqp, int sign)
4672 {
4673         int i;
4674         struct isl_to_poly_data data;
4675
4676         if (sign == 0)
4677                 return pwqp_drop_floors(pwqp);
4678
4679         if (!pwqp)
4680                 return NULL;
4681
4682         data.sign = sign;
4683         data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4684
4685         for (i = 0; i < pwqp->n; ++i) {
4686                 if (pwqp->p[i].qp->div->n_row == 0) {
4687                         isl_pw_qpolynomial *t;
4688                         t = isl_pw_qpolynomial_alloc(
4689                                         isl_set_copy(pwqp->p[i].set),
4690                                         isl_qpolynomial_copy(pwqp->p[i].qp));
4691                         data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4692                         continue;
4693                 }
4694                 data.qp = pwqp->p[i].qp;
4695                 if (isl_set_foreach_orthant(pwqp->p[i].set,
4696                                         &to_polynomial_on_orthant, &data) < 0)
4697                         goto error;
4698         }
4699
4700         isl_pw_qpolynomial_free(pwqp);
4701
4702         return data.res;
4703 error:
4704         isl_pw_qpolynomial_free(pwqp);
4705         isl_pw_qpolynomial_free(data.res);
4706         return NULL;
4707 }
4708
4709 static int poly_entry(void **entry, void *user)
4710 {
4711         int *sign = user;
4712         isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4713
4714         *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4715
4716         return *pwqp ? 0 : -1;
4717 }
4718
4719 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4720         __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4721 {
4722         upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4723         if (!upwqp)
4724                 return NULL;
4725
4726         if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4727                                    &poly_entry, &sign) < 0)
4728                 goto error;
4729
4730         return upwqp;
4731 error:
4732         isl_union_pw_qpolynomial_free(upwqp);
4733         return NULL;
4734 }
4735
4736 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4737         __isl_take isl_qpolynomial *qp)
4738 {
4739         int i, k;
4740         isl_space *dim;
4741         isl_vec *aff = NULL;
4742         isl_basic_map *bmap = NULL;
4743         unsigned pos;
4744         unsigned n_div;
4745
4746         if (!qp)
4747                 return NULL;
4748         if (!isl_upoly_is_affine(qp->upoly))
4749                 isl_die(qp->dim->ctx, isl_error_invalid,
4750                         "input quasi-polynomial not affine", goto error);
4751         aff = isl_qpolynomial_extract_affine(qp);
4752         if (!aff)
4753                 goto error;
4754         dim = isl_qpolynomial_get_space(qp);
4755         pos = 1 + isl_space_offset(dim, isl_dim_out);
4756         n_div = qp->div->n_row;
4757         bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4758
4759         for (i = 0; i < n_div; ++i) {
4760                 k = isl_basic_map_alloc_div(bmap);
4761                 if (k < 0)
4762                         goto error;
4763                 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4764                 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4765                 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4766                         goto error;
4767         }
4768         k = isl_basic_map_alloc_equality(bmap);
4769         if (k < 0)
4770                 goto error;
4771         isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4772         isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4773         isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4774
4775         isl_vec_free(aff);
4776         isl_qpolynomial_free(qp);
4777         bmap = isl_basic_map_finalize(bmap);
4778         return bmap;
4779 error:
4780         isl_vec_free(aff);
4781         isl_qpolynomial_free(qp);
4782         isl_basic_map_free(bmap);
4783         return NULL;
4784 }