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 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2187  */
2188 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2189         __isl_take isl_space *domain, __isl_take isl_val *val)
2190 {
2191         isl_qpolynomial *qp;
2192         struct isl_upoly_cst *cst;
2193
2194         if (!domain || !val)
2195                 goto error;
2196
2197         qp = isl_qpolynomial_alloc(domain, 0, isl_upoly_zero(domain->ctx));
2198         if (!qp)
2199                 goto error;
2200
2201         cst = isl_upoly_as_cst(qp->upoly);
2202         isl_int_set(cst->n, val->n);
2203         isl_int_set(cst->d, val->d);
2204
2205         isl_val_free(val);
2206         return qp;
2207 error:
2208         isl_space_free(domain);
2209         isl_val_free(val);
2210         return NULL;
2211 }
2212
2213 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2214 {
2215         struct isl_upoly_rec *rec;
2216         int i;
2217
2218         if (!up)
2219                 return -1;
2220
2221         if (isl_upoly_is_cst(up))
2222                 return 0;
2223
2224         if (up->var < d)
2225                 active[up->var] = 1;
2226
2227         rec = isl_upoly_as_rec(up);
2228         for (i = 0; i < rec->n; ++i)
2229                 if (up_set_active(rec->p[i], active, d) < 0)
2230                         return -1;
2231
2232         return 0;
2233 }
2234
2235 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2236 {
2237         int i, j;
2238         int d = isl_space_dim(qp->dim, isl_dim_all);
2239
2240         if (!qp || !active)
2241                 return -1;
2242
2243         for (i = 0; i < d; ++i)
2244                 for (j = 0; j < qp->div->n_row; ++j) {
2245                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
2246                                 continue;
2247                         active[i] = 1;
2248                         break;
2249                 }
2250
2251         return up_set_active(qp->upoly, active, d);
2252 }
2253
2254 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2255         enum isl_dim_type type, unsigned first, unsigned n)
2256 {
2257         int i;
2258         int *active = NULL;
2259         int involves = 0;
2260
2261         if (!qp)
2262                 return -1;
2263         if (n == 0)
2264                 return 0;
2265
2266         isl_assert(qp->dim->ctx,
2267                     first + n <= isl_qpolynomial_dim(qp, type), return -1);
2268         isl_assert(qp->dim->ctx, type == isl_dim_param ||
2269                                  type == isl_dim_in, return -1);
2270
2271         active = isl_calloc_array(qp->dim->ctx, int,
2272                                         isl_space_dim(qp->dim, isl_dim_all));
2273         if (set_active(qp, active) < 0)
2274                 goto error;
2275
2276         if (type == isl_dim_in)
2277                 first += isl_space_dim(qp->dim, isl_dim_param);
2278         for (i = 0; i < n; ++i)
2279                 if (active[first + i]) {
2280                         involves = 1;
2281                         break;
2282                 }
2283
2284         free(active);
2285
2286         return involves;
2287 error:
2288         free(active);
2289         return -1;
2290 }
2291
2292 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2293  * of the divs that do appear in the quasi-polynomial.
2294  */
2295 static __isl_give isl_qpolynomial *remove_redundant_divs(
2296         __isl_take isl_qpolynomial *qp)
2297 {
2298         int i, j;
2299         int d;
2300         int len;
2301         int skip;
2302         int *active = NULL;
2303         int *reordering = NULL;
2304         int redundant = 0;
2305         int n_div;
2306         isl_ctx *ctx;
2307
2308         if (!qp)
2309                 return NULL;
2310         if (qp->div->n_row == 0)
2311                 return qp;
2312
2313         d = isl_space_dim(qp->dim, isl_dim_all);
2314         len = qp->div->n_col - 2;
2315         ctx = isl_qpolynomial_get_ctx(qp);
2316         active = isl_calloc_array(ctx, int, len);
2317         if (!active)
2318                 goto error;
2319
2320         if (up_set_active(qp->upoly, active, len) < 0)
2321                 goto error;
2322
2323         for (i = qp->div->n_row - 1; i >= 0; --i) {
2324                 if (!active[d + i]) {
2325                         redundant = 1;
2326                         continue;
2327                 }
2328                 for (j = 0; j < i; ++j) {
2329                         if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2330                                 continue;
2331                         active[d + j] = 1;
2332                         break;
2333                 }
2334         }
2335
2336         if (!redundant) {
2337                 free(active);
2338                 return qp;
2339         }
2340
2341         reordering = isl_alloc_array(qp->div->ctx, int, len);
2342         if (!reordering)
2343                 goto error;
2344
2345         for (i = 0; i < d; ++i)
2346                 reordering[i] = i;
2347
2348         skip = 0;
2349         n_div = qp->div->n_row;
2350         for (i = 0; i < n_div; ++i) {
2351                 if (!active[d + i]) {
2352                         qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2353                         qp->div = isl_mat_drop_cols(qp->div,
2354                                                     2 + d + i - skip, 1);
2355                         skip++;
2356                 }
2357                 reordering[d + i] = d + i - skip;
2358         }
2359
2360         qp->upoly = reorder(qp->upoly, reordering);
2361
2362         if (!qp->upoly || !qp->div)
2363                 goto error;
2364
2365         free(active);
2366         free(reordering);
2367
2368         return qp;
2369 error:
2370         free(active);
2371         free(reordering);
2372         isl_qpolynomial_free(qp);
2373         return NULL;
2374 }
2375
2376 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2377         unsigned first, unsigned n)
2378 {
2379         int i;
2380         struct isl_upoly_rec *rec;
2381
2382         if (!up)
2383                 return NULL;
2384         if (n == 0 || up->var < 0 || up->var < first)
2385                 return up;
2386         if (up->var < first + n) {
2387                 up = replace_by_constant_term(up);
2388                 return isl_upoly_drop(up, first, n);
2389         }
2390         up = isl_upoly_cow(up);
2391         if (!up)
2392                 return NULL;
2393         up->var -= n;
2394         rec = isl_upoly_as_rec(up);
2395         if (!rec)
2396                 goto error;
2397
2398         for (i = 0; i < rec->n; ++i) {
2399                 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2400                 if (!rec->p[i])
2401                         goto error;
2402         }
2403
2404         return up;
2405 error:
2406         isl_upoly_free(up);
2407         return NULL;
2408 }
2409
2410 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2411         __isl_take isl_qpolynomial *qp,
2412         enum isl_dim_type type, unsigned pos, const char *s)
2413 {
2414         qp = isl_qpolynomial_cow(qp);
2415         if (!qp)
2416                 return NULL;
2417         qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2418         if (!qp->dim)
2419                 goto error;
2420         return qp;
2421 error:
2422         isl_qpolynomial_free(qp);
2423         return NULL;
2424 }
2425
2426 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2427         __isl_take isl_qpolynomial *qp,
2428         enum isl_dim_type type, unsigned first, unsigned n)
2429 {
2430         if (!qp)
2431                 return NULL;
2432         if (type == isl_dim_out)
2433                 isl_die(qp->dim->ctx, isl_error_invalid,
2434                         "cannot drop output/set dimension",
2435                         goto error);
2436         if (type == isl_dim_in)
2437                 type = isl_dim_set;
2438         if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2439                 return qp;
2440
2441         qp = isl_qpolynomial_cow(qp);
2442         if (!qp)
2443                 return NULL;
2444
2445         isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2446                         goto error);
2447         isl_assert(qp->dim->ctx, type == isl_dim_param ||
2448                                  type == isl_dim_set, goto error);
2449
2450         qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2451         if (!qp->dim)
2452                 goto error;
2453
2454         if (type == isl_dim_set)
2455                 first += isl_space_dim(qp->dim, isl_dim_param);
2456
2457         qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2458         if (!qp->div)
2459                 goto error;
2460
2461         qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2462         if (!qp->upoly)
2463                 goto error;
2464
2465         return qp;
2466 error:
2467         isl_qpolynomial_free(qp);
2468         return NULL;
2469 }
2470
2471 /* Project the domain of the quasi-polynomial onto its parameter space.
2472  * The quasi-polynomial may not involve any of the domain dimensions.
2473  */
2474 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2475         __isl_take isl_qpolynomial *qp)
2476 {
2477         isl_space *space;
2478         unsigned n;
2479         int involves;
2480
2481         n = isl_qpolynomial_dim(qp, isl_dim_in);
2482         involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2483         if (involves < 0)
2484                 return isl_qpolynomial_free(qp);
2485         if (involves)
2486                 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2487                         "polynomial involves some of the domain dimensions",
2488                         return isl_qpolynomial_free(qp));
2489         qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2490         space = isl_qpolynomial_get_domain_space(qp);
2491         space = isl_space_params(space);
2492         qp = isl_qpolynomial_reset_domain_space(qp, space);
2493         return qp;
2494 }
2495
2496 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2497         __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2498 {
2499         int i, j, k;
2500         isl_int denom;
2501         unsigned total;
2502         unsigned n_div;
2503         struct isl_upoly *up;
2504
2505         if (!eq)
2506                 goto error;
2507         if (eq->n_eq == 0) {
2508                 isl_basic_set_free(eq);
2509                 return qp;
2510         }
2511
2512         qp = isl_qpolynomial_cow(qp);
2513         if (!qp)
2514                 goto error;
2515         qp->div = isl_mat_cow(qp->div);
2516         if (!qp->div)
2517                 goto error;
2518
2519         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2520         n_div = eq->n_div;
2521         isl_int_init(denom);
2522         for (i = 0; i < eq->n_eq; ++i) {
2523                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2524                 if (j < 0 || j == 0 || j >= total)
2525                         continue;
2526
2527                 for (k = 0; k < qp->div->n_row; ++k) {
2528                         if (isl_int_is_zero(qp->div->row[k][1 + j]))
2529                                 continue;
2530                         isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2531                                         &qp->div->row[k][0]);
2532                         normalize_div(qp, k);
2533                 }
2534
2535                 if (isl_int_is_pos(eq->eq[i][j]))
2536                         isl_seq_neg(eq->eq[i], eq->eq[i], total);
2537                 isl_int_abs(denom, eq->eq[i][j]);
2538                 isl_int_set_si(eq->eq[i][j], 0);
2539
2540                 up = isl_upoly_from_affine(qp->dim->ctx,
2541                                                    eq->eq[i], denom, total);
2542                 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2543                 isl_upoly_free(up);
2544         }
2545         isl_int_clear(denom);
2546
2547         if (!qp->upoly)
2548                 goto error;
2549
2550         isl_basic_set_free(eq);
2551
2552         qp = substitute_non_divs(qp);
2553         qp = sort_divs(qp);
2554
2555         return qp;
2556 error:
2557         isl_basic_set_free(eq);
2558         isl_qpolynomial_free(qp);
2559         return NULL;
2560 }
2561
2562 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2563  */
2564 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2565         __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2566 {
2567         if (!qp || !eq)
2568                 goto error;
2569         if (qp->div->n_row > 0)
2570                 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2571         return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2572 error:
2573         isl_basic_set_free(eq);
2574         isl_qpolynomial_free(qp);
2575         return NULL;
2576 }
2577
2578 static __isl_give isl_basic_set *add_div_constraints(
2579         __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2580 {
2581         int i;
2582         unsigned total;
2583
2584         if (!bset || !div)
2585                 goto error;
2586
2587         bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2588         if (!bset)
2589                 goto error;
2590         total = isl_basic_set_total_dim(bset);
2591         for (i = 0; i < div->n_row; ++i)
2592                 if (isl_basic_set_add_div_constraints_var(bset,
2593                                     total - div->n_row + i, div->row[i]) < 0)
2594                         goto error;
2595
2596         isl_mat_free(div);
2597         return bset;
2598 error:
2599         isl_mat_free(div);
2600         isl_basic_set_free(bset);
2601         return NULL;
2602 }
2603
2604 /* Look for equalities among the variables shared by context and qp
2605  * and the integer divisions of qp, if any.
2606  * The equalities are then used to eliminate variables and/or integer
2607  * divisions from qp.
2608  */
2609 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2610         __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2611 {
2612         isl_basic_set *aff;
2613
2614         if (!qp)
2615                 goto error;
2616         if (qp->div->n_row > 0) {
2617                 isl_basic_set *bset;
2618                 context = isl_set_add_dims(context, isl_dim_set,
2619                                             qp->div->n_row);
2620                 bset = isl_basic_set_universe(isl_set_get_space(context));
2621                 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2622                 context = isl_set_intersect(context,
2623                                             isl_set_from_basic_set(bset));
2624         }
2625
2626         aff = isl_set_affine_hull(context);
2627         return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2628 error:
2629         isl_qpolynomial_free(qp);
2630         isl_set_free(context);
2631         return NULL;
2632 }
2633
2634 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2635         __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2636 {
2637         isl_space *space = isl_qpolynomial_get_domain_space(qp);
2638         isl_set *dom_context = isl_set_universe(space);
2639         dom_context = isl_set_intersect_params(dom_context, context);
2640         return isl_qpolynomial_gist(qp, dom_context);
2641 }
2642
2643 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2644         __isl_take isl_qpolynomial *qp)
2645 {
2646         isl_set *dom;
2647
2648         if (!qp)
2649                 return NULL;
2650         if (isl_qpolynomial_is_zero(qp)) {
2651                 isl_space *dim = isl_qpolynomial_get_space(qp);
2652                 isl_qpolynomial_free(qp);
2653                 return isl_pw_qpolynomial_zero(dim);
2654         }
2655
2656         dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2657         return isl_pw_qpolynomial_alloc(dom, qp);
2658 }
2659
2660 #undef PW
2661 #define PW isl_pw_qpolynomial
2662 #undef EL
2663 #define EL isl_qpolynomial
2664 #undef EL_IS_ZERO
2665 #define EL_IS_ZERO is_zero
2666 #undef ZERO
2667 #define ZERO zero
2668 #undef IS_ZERO
2669 #define IS_ZERO is_zero
2670 #undef FIELD
2671 #define FIELD qp
2672 #undef DEFAULT_IS_ZERO
2673 #define DEFAULT_IS_ZERO 1
2674
2675 #define NO_PULLBACK
2676
2677 #include <isl_pw_templ.c>
2678
2679 #undef UNION
2680 #define UNION isl_union_pw_qpolynomial
2681 #undef PART
2682 #define PART isl_pw_qpolynomial
2683 #undef PARTS
2684 #define PARTS pw_qpolynomial
2685 #define ALIGN_DOMAIN
2686
2687 #include <isl_union_templ.c>
2688
2689 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2690 {
2691         if (!pwqp)
2692                 return -1;
2693
2694         if (pwqp->n != -1)
2695                 return 0;
2696
2697         if (!isl_set_plain_is_universe(pwqp->p[0].set))
2698                 return 0;
2699
2700         return isl_qpolynomial_is_one(pwqp->p[0].qp);
2701 }
2702
2703 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2704         __isl_take isl_pw_qpolynomial *pwqp1,
2705         __isl_take isl_pw_qpolynomial *pwqp2)
2706 {
2707         return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2708 }
2709
2710 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2711         __isl_take isl_pw_qpolynomial *pwqp1,
2712         __isl_take isl_pw_qpolynomial *pwqp2)
2713 {
2714         int i, j, n;
2715         struct isl_pw_qpolynomial *res;
2716
2717         if (!pwqp1 || !pwqp2)
2718                 goto error;
2719
2720         isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2721                         goto error);
2722
2723         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2724                 isl_pw_qpolynomial_free(pwqp2);
2725                 return pwqp1;
2726         }
2727
2728         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2729                 isl_pw_qpolynomial_free(pwqp1);
2730                 return pwqp2;
2731         }
2732
2733         if (isl_pw_qpolynomial_is_one(pwqp1)) {
2734                 isl_pw_qpolynomial_free(pwqp1);
2735                 return pwqp2;
2736         }
2737
2738         if (isl_pw_qpolynomial_is_one(pwqp2)) {
2739                 isl_pw_qpolynomial_free(pwqp2);
2740                 return pwqp1;
2741         }
2742
2743         n = pwqp1->n * pwqp2->n;
2744         res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2745
2746         for (i = 0; i < pwqp1->n; ++i) {
2747                 for (j = 0; j < pwqp2->n; ++j) {
2748                         struct isl_set *common;
2749                         struct isl_qpolynomial *prod;
2750                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2751                                                 isl_set_copy(pwqp2->p[j].set));
2752                         if (isl_set_plain_is_empty(common)) {
2753                                 isl_set_free(common);
2754                                 continue;
2755                         }
2756
2757                         prod = isl_qpolynomial_mul(
2758                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
2759                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
2760
2761                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
2762                 }
2763         }
2764
2765         isl_pw_qpolynomial_free(pwqp1);
2766         isl_pw_qpolynomial_free(pwqp2);
2767
2768         return res;
2769 error:
2770         isl_pw_qpolynomial_free(pwqp1);
2771         isl_pw_qpolynomial_free(pwqp2);
2772         return NULL;
2773 }
2774
2775 __isl_give struct isl_upoly *isl_upoly_eval(
2776         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
2777 {
2778         int i;
2779         struct isl_upoly_rec *rec;
2780         struct isl_upoly *res;
2781         struct isl_upoly *base;
2782
2783         if (isl_upoly_is_cst(up)) {
2784                 isl_vec_free(vec);
2785                 return up;
2786         }
2787
2788         rec = isl_upoly_as_rec(up);
2789         if (!rec)
2790                 goto error;
2791
2792         isl_assert(up->ctx, rec->n >= 1, goto error);
2793
2794         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
2795
2796         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2797                                 isl_vec_copy(vec));
2798
2799         for (i = rec->n - 2; i >= 0; --i) {
2800                 res = isl_upoly_mul(res, isl_upoly_copy(base));
2801                 res = isl_upoly_sum(res,
2802                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2803                                                             isl_vec_copy(vec)));
2804         }
2805
2806         isl_upoly_free(base);
2807         isl_upoly_free(up);
2808         isl_vec_free(vec);
2809         return res;
2810 error:
2811         isl_upoly_free(up);
2812         isl_vec_free(vec);
2813         return NULL;
2814 }
2815
2816 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
2817         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
2818 {
2819         isl_vec *ext;
2820         struct isl_upoly *up;
2821         isl_space *dim;
2822
2823         if (!qp || !pnt)
2824                 goto error;
2825         isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2826
2827         if (qp->div->n_row == 0)
2828                 ext = isl_vec_copy(pnt->vec);
2829         else {
2830                 int i;
2831                 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2832                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2833                 if (!ext)
2834                         goto error;
2835
2836                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2837                 for (i = 0; i < qp->div->n_row; ++i) {
2838                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2839                                                 1 + dim + i, &ext->el[1+dim+i]);
2840                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2841                                         qp->div->row[i][0]);
2842                 }
2843         }
2844
2845         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2846         if (!up)
2847                 goto error;
2848
2849         dim = isl_space_copy(qp->dim);
2850         isl_qpolynomial_free(qp);
2851         isl_point_free(pnt);
2852
2853         return isl_qpolynomial_alloc(dim, 0, up);
2854 error:
2855         isl_qpolynomial_free(qp);
2856         isl_point_free(pnt);
2857         return NULL;
2858 }
2859
2860 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2861         __isl_keep struct isl_upoly_cst *cst2)
2862 {
2863         int cmp;
2864         isl_int t;
2865         isl_int_init(t);
2866         isl_int_mul(t, cst1->n, cst2->d);
2867         isl_int_submul(t, cst2->n, cst1->d);
2868         cmp = isl_int_sgn(t);
2869         isl_int_clear(t);
2870         return cmp;
2871 }
2872
2873 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial *qp1,
2874         __isl_keep isl_qpolynomial *qp2)
2875 {
2876         struct isl_upoly_cst *cst1, *cst2;
2877
2878         if (!qp1 || !qp2)
2879                 return -1;
2880         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), return -1);
2881         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), return -1);
2882         if (isl_qpolynomial_is_nan(qp1))
2883                 return -1;
2884         if (isl_qpolynomial_is_nan(qp2))
2885                 return -1;
2886         cst1 = isl_upoly_as_cst(qp1->upoly);
2887         cst2 = isl_upoly_as_cst(qp2->upoly);
2888
2889         return isl_upoly_cmp(cst1, cst2) <= 0;
2890 }
2891
2892 __isl_give isl_qpolynomial *isl_qpolynomial_min_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_max_cst(
2920         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2921 {
2922         struct isl_upoly_cst *cst1, *cst2;
2923         int cmp;
2924
2925         if (!qp1 || !qp2)
2926                 goto error;
2927         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2928         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2929         cst1 = isl_upoly_as_cst(qp1->upoly);
2930         cst2 = isl_upoly_as_cst(qp2->upoly);
2931         cmp = isl_upoly_cmp(cst1, cst2);
2932
2933         if (cmp >= 0) {
2934                 isl_qpolynomial_free(qp2);
2935         } else {
2936                 isl_qpolynomial_free(qp1);
2937                 qp1 = qp2;
2938         }
2939         return qp1;
2940 error:
2941         isl_qpolynomial_free(qp1);
2942         isl_qpolynomial_free(qp2);
2943         return NULL;
2944 }
2945
2946 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2947         __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2948         unsigned first, unsigned n)
2949 {
2950         unsigned total;
2951         unsigned g_pos;
2952         int *exp;
2953
2954         if (!qp)
2955                 return NULL;
2956         if (type == isl_dim_out)
2957                 isl_die(qp->div->ctx, isl_error_invalid,
2958                         "cannot insert output/set dimensions",
2959                         goto error);
2960         if (type == isl_dim_in)
2961                 type = isl_dim_set;
2962         if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2963                 return qp;
2964
2965         qp = isl_qpolynomial_cow(qp);
2966         if (!qp)
2967                 return NULL;
2968
2969         isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
2970                     goto error);
2971
2972         g_pos = pos(qp->dim, type) + first;
2973
2974         qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
2975         if (!qp->div)
2976                 goto error;
2977
2978         total = qp->div->n_col - 2;
2979         if (total > g_pos) {
2980                 int i;
2981                 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
2982                 if (!exp)
2983                         goto error;
2984                 for (i = 0; i < total - g_pos; ++i)
2985                         exp[i] = i + n;
2986                 qp->upoly = expand(qp->upoly, exp, g_pos);
2987                 free(exp);
2988                 if (!qp->upoly)
2989                         goto error;
2990         }
2991
2992         qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
2993         if (!qp->dim)
2994                 goto error;
2995
2996         return qp;
2997 error:
2998         isl_qpolynomial_free(qp);
2999         return NULL;
3000 }
3001
3002 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3003         __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3004 {
3005         unsigned pos;
3006
3007         pos = isl_qpolynomial_dim(qp, type);
3008
3009         return isl_qpolynomial_insert_dims(qp, type, pos, n);
3010 }
3011
3012 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3013         __isl_take isl_pw_qpolynomial *pwqp,
3014         enum isl_dim_type type, unsigned n)
3015 {
3016         unsigned pos;
3017
3018         pos = isl_pw_qpolynomial_dim(pwqp, type);
3019
3020         return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3021 }
3022
3023 static int *reordering_move(isl_ctx *ctx,
3024         unsigned len, unsigned dst, unsigned src, unsigned n)
3025 {
3026         int i;
3027         int *reordering;
3028
3029         reordering = isl_alloc_array(ctx, int, len);
3030         if (!reordering)
3031                 return NULL;
3032
3033         if (dst <= src) {
3034                 for (i = 0; i < dst; ++i)
3035                         reordering[i] = i;
3036                 for (i = 0; i < n; ++i)
3037                         reordering[src + i] = dst + i;
3038                 for (i = 0; i < src - dst; ++i)
3039                         reordering[dst + i] = dst + n + i;
3040                 for (i = 0; i < len - src - n; ++i)
3041                         reordering[src + n + i] = src + n + i;
3042         } else {
3043                 for (i = 0; i < src; ++i)
3044                         reordering[i] = i;
3045                 for (i = 0; i < n; ++i)
3046                         reordering[src + i] = dst + i;
3047                 for (i = 0; i < dst - src; ++i)
3048                         reordering[src + n + i] = src + i;
3049                 for (i = 0; i < len - dst - n; ++i)
3050                         reordering[dst + n + i] = dst + n + i;
3051         }
3052
3053         return reordering;
3054 }
3055
3056 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3057         __isl_take isl_qpolynomial *qp,
3058         enum isl_dim_type dst_type, unsigned dst_pos,
3059         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3060 {
3061         unsigned g_dst_pos;
3062         unsigned g_src_pos;
3063         int *reordering;
3064
3065         qp = isl_qpolynomial_cow(qp);
3066         if (!qp)
3067                 return NULL;
3068
3069         if (dst_type == isl_dim_out || src_type == isl_dim_out)
3070                 isl_die(qp->dim->ctx, isl_error_invalid,
3071                         "cannot move output/set dimension",
3072                         goto error);
3073         if (dst_type == isl_dim_in)
3074                 dst_type = isl_dim_set;
3075         if (src_type == isl_dim_in)
3076                 src_type = isl_dim_set;
3077
3078         isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3079                 goto error);
3080
3081         g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3082         g_src_pos = pos(qp->dim, src_type) + src_pos;
3083         if (dst_type > src_type)
3084                 g_dst_pos -= n;
3085
3086         qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3087         if (!qp->div)
3088                 goto error;
3089         qp = sort_divs(qp);
3090         if (!qp)
3091                 goto error;
3092
3093         reordering = reordering_move(qp->dim->ctx,
3094                                 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3095         if (!reordering)
3096                 goto error;
3097
3098         qp->upoly = reorder(qp->upoly, reordering);
3099         free(reordering);
3100         if (!qp->upoly)
3101                 goto error;
3102
3103         qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3104         if (!qp->dim)
3105                 goto error;
3106
3107         return qp;
3108 error:
3109         isl_qpolynomial_free(qp);
3110         return NULL;
3111 }
3112
3113 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3114         isl_int *f, isl_int denom)
3115 {
3116         struct isl_upoly *up;
3117
3118         dim = isl_space_domain(dim);
3119         if (!dim)
3120                 return NULL;
3121
3122         up = isl_upoly_from_affine(dim->ctx, f, denom,
3123                                         1 + isl_space_dim(dim, isl_dim_all));
3124
3125         return isl_qpolynomial_alloc(dim, 0, up);
3126 }
3127
3128 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3129 {
3130         isl_ctx *ctx;
3131         struct isl_upoly *up;
3132         isl_qpolynomial *qp;
3133
3134         if (!aff)
3135                 return NULL;
3136
3137         ctx = isl_aff_get_ctx(aff);
3138         up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3139                                     aff->v->size - 1);
3140
3141         qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3142                                     aff->ls->div->n_row, up);
3143         if (!qp)
3144                 goto error;
3145
3146         isl_mat_free(qp->div);
3147         qp->div = isl_mat_copy(aff->ls->div);
3148         qp->div = isl_mat_cow(qp->div);
3149         if (!qp->div)
3150                 goto error;
3151
3152         isl_aff_free(aff);
3153         qp = reduce_divs(qp);
3154         qp = remove_redundant_divs(qp);
3155         return qp;
3156 error:
3157         isl_aff_free(aff);
3158         return NULL;
3159 }
3160
3161 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3162         __isl_take isl_pw_aff *pwaff)
3163 {
3164         int i;
3165         isl_pw_qpolynomial *pwqp;
3166
3167         if (!pwaff)
3168                 return NULL;
3169
3170         pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3171                                                 pwaff->n);
3172
3173         for (i = 0; i < pwaff->n; ++i) {
3174                 isl_set *dom;
3175                 isl_qpolynomial *qp;
3176
3177                 dom = isl_set_copy(pwaff->p[i].set);
3178                 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3179                 pwqp = isl_pw_qpolynomial_add_piece(pwqp,  dom, qp);
3180         }
3181
3182         isl_pw_aff_free(pwaff);
3183         return pwqp;
3184 }
3185
3186 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3187         __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3188 {
3189         isl_aff *aff;
3190
3191         aff = isl_constraint_get_bound(c, type, pos);
3192         isl_constraint_free(c);
3193         return isl_qpolynomial_from_aff(aff);
3194 }
3195
3196 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3197  * in "qp" by subs[i].
3198  */
3199 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3200         __isl_take isl_qpolynomial *qp,
3201         enum isl_dim_type type, unsigned first, unsigned n,
3202         __isl_keep isl_qpolynomial **subs)
3203 {
3204         int i;
3205         struct isl_upoly **ups;
3206
3207         if (n == 0)
3208                 return qp;
3209
3210         qp = isl_qpolynomial_cow(qp);
3211         if (!qp)
3212                 return NULL;
3213
3214         if (type == isl_dim_out)
3215                 isl_die(qp->dim->ctx, isl_error_invalid,
3216                         "cannot substitute output/set dimension",
3217                         goto error);
3218         if (type == isl_dim_in)
3219                 type = isl_dim_set;
3220
3221         for (i = 0; i < n; ++i)
3222                 if (!subs[i])
3223                         goto error;
3224
3225         isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3226                         goto error);
3227
3228         for (i = 0; i < n; ++i)
3229                 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3230                                 goto error);
3231
3232         isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3233         for (i = 0; i < n; ++i)
3234                 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3235
3236         first += pos(qp->dim, type);
3237
3238         ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3239         if (!ups)
3240                 goto error;
3241         for (i = 0; i < n; ++i)
3242                 ups[i] = subs[i]->upoly;
3243
3244         qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3245
3246         free(ups);
3247
3248         if (!qp->upoly)
3249                 goto error;
3250
3251         return qp;
3252 error:
3253         isl_qpolynomial_free(qp);
3254         return NULL;
3255 }
3256
3257 /* Extend "bset" with extra set dimensions for each integer division
3258  * in "qp" and then call "fn" with the extended bset and the polynomial
3259  * that results from replacing each of the integer divisions by the
3260  * corresponding extra set dimension.
3261  */
3262 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3263         __isl_keep isl_basic_set *bset,
3264         int (*fn)(__isl_take isl_basic_set *bset,
3265                   __isl_take isl_qpolynomial *poly, void *user), void *user)
3266 {
3267         isl_space *dim;
3268         isl_mat *div;
3269         isl_qpolynomial *poly;
3270
3271         if (!qp || !bset)
3272                 goto error;
3273         if (qp->div->n_row == 0)
3274                 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3275                           user);
3276
3277         div = isl_mat_copy(qp->div);
3278         dim = isl_space_copy(qp->dim);
3279         dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3280         poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3281         bset = isl_basic_set_copy(bset);
3282         bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3283         bset = add_div_constraints(bset, div);
3284
3285         return fn(bset, poly, user);
3286 error:
3287         return -1;
3288 }
3289
3290 /* Return total degree in variables first (inclusive) up to last (exclusive).
3291  */
3292 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3293 {
3294         int deg = -1;
3295         int i;
3296         struct isl_upoly_rec *rec;
3297
3298         if (!up)
3299                 return -2;
3300         if (isl_upoly_is_zero(up))
3301                 return -1;
3302         if (isl_upoly_is_cst(up) || up->var < first)
3303                 return 0;
3304
3305         rec = isl_upoly_as_rec(up);
3306         if (!rec)
3307                 return -2;
3308
3309         for (i = 0; i < rec->n; ++i) {
3310                 int d;
3311
3312                 if (isl_upoly_is_zero(rec->p[i]))
3313                         continue;
3314                 d = isl_upoly_degree(rec->p[i], first, last);
3315                 if (up->var < last)
3316                         d += i;
3317                 if (d > deg)
3318                         deg = d;
3319         }
3320
3321         return deg;
3322 }
3323
3324 /* Return total degree in set variables.
3325  */
3326 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3327 {
3328         unsigned ovar;
3329         unsigned nvar;
3330
3331         if (!poly)
3332                 return -2;
3333
3334         ovar = isl_space_offset(poly->dim, isl_dim_set);
3335         nvar = isl_space_dim(poly->dim, isl_dim_set);
3336         return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3337 }
3338
3339 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3340         unsigned pos, int deg)
3341 {
3342         int i;
3343         struct isl_upoly_rec *rec;
3344
3345         if (!up)
3346                 return NULL;
3347
3348         if (isl_upoly_is_cst(up) || up->var < pos) {
3349                 if (deg == 0)
3350                         return isl_upoly_copy(up);
3351                 else
3352                         return isl_upoly_zero(up->ctx);
3353         }
3354
3355         rec = isl_upoly_as_rec(up);
3356         if (!rec)
3357                 return NULL;
3358
3359         if (up->var == pos) {
3360                 if (deg < rec->n)
3361                         return isl_upoly_copy(rec->p[deg]);
3362                 else
3363                         return isl_upoly_zero(up->ctx);
3364         }
3365
3366         up = isl_upoly_copy(up);
3367         up = isl_upoly_cow(up);
3368         rec = isl_upoly_as_rec(up);
3369         if (!rec)
3370                 goto error;
3371
3372         for (i = 0; i < rec->n; ++i) {
3373                 struct isl_upoly *t;
3374                 t = isl_upoly_coeff(rec->p[i], pos, deg);
3375                 if (!t)
3376                         goto error;
3377                 isl_upoly_free(rec->p[i]);
3378                 rec->p[i] = t;
3379         }
3380
3381         return up;
3382 error:
3383         isl_upoly_free(up);
3384         return NULL;
3385 }
3386
3387 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3388  */
3389 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3390         __isl_keep isl_qpolynomial *qp,
3391         enum isl_dim_type type, unsigned t_pos, int deg)
3392 {
3393         unsigned g_pos;
3394         struct isl_upoly *up;
3395         isl_qpolynomial *c;
3396
3397         if (!qp)
3398                 return NULL;
3399
3400         if (type == isl_dim_out)
3401                 isl_die(qp->div->ctx, isl_error_invalid,
3402                         "output/set dimension does not have a coefficient",
3403                         return NULL);
3404         if (type == isl_dim_in)
3405                 type = isl_dim_set;
3406
3407         isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3408                         return NULL);
3409
3410         g_pos = pos(qp->dim, type) + t_pos;
3411         up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3412
3413         c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3414         if (!c)
3415                 return NULL;
3416         isl_mat_free(c->div);
3417         c->div = isl_mat_copy(qp->div);
3418         if (!c->div)
3419                 goto error;
3420         return c;
3421 error:
3422         isl_qpolynomial_free(c);
3423         return NULL;
3424 }
3425
3426 /* Homogenize the polynomial in the variables first (inclusive) up to
3427  * last (exclusive) by inserting powers of variable first.
3428  * Variable first is assumed not to appear in the input.
3429  */
3430 __isl_give struct isl_upoly *isl_upoly_homogenize(
3431         __isl_take struct isl_upoly *up, int deg, int target,
3432         int first, int last)
3433 {
3434         int i;
3435         struct isl_upoly_rec *rec;
3436
3437         if (!up)
3438                 return NULL;
3439         if (isl_upoly_is_zero(up))
3440                 return up;
3441         if (deg == target)
3442                 return up;
3443         if (isl_upoly_is_cst(up) || up->var < first) {
3444                 struct isl_upoly *hom;
3445
3446                 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3447                 if (!hom)
3448                         goto error;
3449                 rec = isl_upoly_as_rec(hom);
3450                 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3451
3452                 return hom;
3453         }
3454
3455         up = isl_upoly_cow(up);
3456         rec = isl_upoly_as_rec(up);
3457         if (!rec)
3458                 goto error;
3459
3460         for (i = 0; i < rec->n; ++i) {
3461                 if (isl_upoly_is_zero(rec->p[i]))
3462                         continue;
3463                 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3464                                 up->var < last ? deg + i : i, target,
3465                                 first, last);
3466                 if (!rec->p[i])
3467                         goto error;
3468         }
3469
3470         return up;
3471 error:
3472         isl_upoly_free(up);
3473         return NULL;
3474 }
3475
3476 /* Homogenize the polynomial in the set variables by introducing
3477  * powers of an extra set variable at position 0.
3478  */
3479 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3480         __isl_take isl_qpolynomial *poly)
3481 {
3482         unsigned ovar;
3483         unsigned nvar;
3484         int deg = isl_qpolynomial_degree(poly);
3485
3486         if (deg < -1)
3487                 goto error;
3488
3489         poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3490         poly = isl_qpolynomial_cow(poly);
3491         if (!poly)
3492                 goto error;
3493
3494         ovar = isl_space_offset(poly->dim, isl_dim_set);
3495         nvar = isl_space_dim(poly->dim, isl_dim_set);
3496         poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3497                                                 ovar, ovar + nvar);
3498         if (!poly->upoly)
3499                 goto error;
3500
3501         return poly;
3502 error:
3503         isl_qpolynomial_free(poly);
3504         return NULL;
3505 }
3506
3507 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3508         __isl_take isl_mat *div)
3509 {
3510         isl_term *term;
3511         int n;
3512
3513         if (!dim || !div)
3514                 goto error;
3515
3516         n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3517
3518         term = isl_calloc(dim->ctx, struct isl_term,
3519                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
3520         if (!term)
3521                 goto error;
3522
3523         term->ref = 1;
3524         term->dim = dim;
3525         term->div = div;
3526         isl_int_init(term->n);
3527         isl_int_init(term->d);
3528
3529         return term;
3530 error:
3531         isl_space_free(dim);
3532         isl_mat_free(div);
3533         return NULL;
3534 }
3535
3536 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3537 {
3538         if (!term)
3539                 return NULL;
3540
3541         term->ref++;
3542         return term;
3543 }
3544
3545 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3546 {
3547         int i;
3548         isl_term *dup;
3549         unsigned total;
3550
3551         if (!term)
3552                 return NULL;
3553
3554         total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3555
3556         dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3557         if (!dup)
3558                 return NULL;
3559
3560         isl_int_set(dup->n, term->n);
3561         isl_int_set(dup->d, term->d);
3562
3563         for (i = 0; i < total; ++i)
3564                 dup->pow[i] = term->pow[i];
3565
3566         return dup;
3567 }
3568
3569 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3570 {
3571         if (!term)
3572                 return NULL;
3573
3574         if (term->ref == 1)
3575                 return term;
3576         term->ref--;
3577         return isl_term_dup(term);
3578 }
3579
3580 void isl_term_free(__isl_take isl_term *term)
3581 {
3582         if (!term)
3583                 return;
3584
3585         if (--term->ref > 0)
3586                 return;
3587
3588         isl_space_free(term->dim);
3589         isl_mat_free(term->div);
3590         isl_int_clear(term->n);
3591         isl_int_clear(term->d);
3592         free(term);
3593 }
3594
3595 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3596 {
3597         if (!term)
3598                 return 0;
3599
3600         switch (type) {
3601         case isl_dim_param:
3602         case isl_dim_in:
3603         case isl_dim_out:       return isl_space_dim(term->dim, type);
3604         case isl_dim_div:       return term->div->n_row;
3605         case isl_dim_all:       return isl_space_dim(term->dim, isl_dim_all) +
3606                                                                 term->div->n_row;
3607         default:                return 0;
3608         }
3609 }
3610
3611 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3612 {
3613         return term ? term->dim->ctx : NULL;
3614 }
3615
3616 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3617 {
3618         if (!term)
3619                 return;
3620         isl_int_set(*n, term->n);
3621 }
3622
3623 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3624 {
3625         if (!term)
3626                 return;
3627         isl_int_set(*d, term->d);
3628 }
3629
3630 /* Return the coefficient of the term "term".
3631  */
3632 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3633 {
3634         if (!term)
3635                 return NULL;
3636
3637         return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3638                                         term->n, term->d);
3639 }
3640
3641 int isl_term_get_exp(__isl_keep isl_term *term,
3642         enum isl_dim_type type, unsigned pos)
3643 {
3644         if (!term)
3645                 return -1;
3646
3647         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3648
3649         if (type >= isl_dim_set)
3650                 pos += isl_space_dim(term->dim, isl_dim_param);
3651         if (type >= isl_dim_div)
3652                 pos += isl_space_dim(term->dim, isl_dim_set);
3653
3654         return term->pow[pos];
3655 }
3656
3657 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3658 {
3659         isl_local_space *ls;
3660         isl_aff *aff;
3661
3662         if (!term)
3663                 return NULL;
3664
3665         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3666                         return NULL);
3667
3668         ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3669                                         isl_mat_copy(term->div));
3670         aff = isl_aff_alloc(ls);
3671         if (!aff)
3672                 return NULL;
3673
3674         isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3675
3676         aff = isl_aff_normalize(aff);
3677
3678         return aff;
3679 }
3680
3681 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3682         int (*fn)(__isl_take isl_term *term, void *user),
3683         __isl_take isl_term *term, void *user)
3684 {
3685         int i;
3686         struct isl_upoly_rec *rec;
3687
3688         if (!up || !term)
3689                 goto error;
3690
3691         if (isl_upoly_is_zero(up))
3692                 return term;
3693
3694         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3695         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3696         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3697
3698         if (isl_upoly_is_cst(up)) {
3699                 struct isl_upoly_cst *cst;
3700                 cst = isl_upoly_as_cst(up);
3701                 if (!cst)
3702                         goto error;
3703                 term = isl_term_cow(term);
3704                 if (!term)
3705                         goto error;
3706                 isl_int_set(term->n, cst->n);
3707                 isl_int_set(term->d, cst->d);
3708                 if (fn(isl_term_copy(term), user) < 0)
3709                         goto error;
3710                 return term;
3711         }
3712
3713         rec = isl_upoly_as_rec(up);
3714         if (!rec)
3715                 goto error;
3716
3717         for (i = 0; i < rec->n; ++i) {
3718                 term = isl_term_cow(term);
3719                 if (!term)
3720                         goto error;
3721                 term->pow[up->var] = i;
3722                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3723                 if (!term)
3724                         goto error;
3725         }
3726         term->pow[up->var] = 0;
3727
3728         return term;
3729 error:
3730         isl_term_free(term);
3731         return NULL;
3732 }
3733
3734 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3735         int (*fn)(__isl_take isl_term *term, void *user), void *user)
3736 {
3737         isl_term *term;
3738
3739         if (!qp)
3740                 return -1;
3741
3742         term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3743         if (!term)
3744                 return -1;
3745
3746         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3747
3748         isl_term_free(term);
3749
3750         return term ? 0 : -1;
3751 }
3752
3753 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3754 {
3755         struct isl_upoly *up;
3756         isl_qpolynomial *qp;
3757         int i, n;
3758
3759         if (!term)
3760                 return NULL;
3761
3762         n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3763
3764         up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3765         for (i = 0; i < n; ++i) {
3766                 if (!term->pow[i])
3767                         continue;
3768                 up = isl_upoly_mul(up,
3769                         isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3770         }
3771
3772         qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3773         if (!qp)
3774                 goto error;
3775         isl_mat_free(qp->div);
3776         qp->div = isl_mat_copy(term->div);
3777         if (!qp->div)
3778                 goto error;
3779
3780         isl_term_free(term);
3781         return qp;
3782 error:
3783         isl_qpolynomial_free(qp);
3784         isl_term_free(term);
3785         return NULL;
3786 }
3787
3788 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3789         __isl_take isl_space *dim)
3790 {
3791         int i;
3792         int extra;
3793         unsigned total;
3794
3795         if (!qp || !dim)
3796                 goto error;
3797
3798         if (isl_space_is_equal(qp->dim, dim)) {
3799                 isl_space_free(dim);
3800                 return qp;
3801         }
3802
3803         qp = isl_qpolynomial_cow(qp);
3804         if (!qp)
3805                 goto error;
3806
3807         extra = isl_space_dim(dim, isl_dim_set) -
3808                         isl_space_dim(qp->dim, isl_dim_set);
3809         total = isl_space_dim(qp->dim, isl_dim_all);
3810         if (qp->div->n_row) {
3811                 int *exp;
3812
3813                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3814                 if (!exp)
3815                         goto error;
3816                 for (i = 0; i < qp->div->n_row; ++i)
3817                         exp[i] = extra + i;
3818                 qp->upoly = expand(qp->upoly, exp, total);
3819                 free(exp);
3820                 if (!qp->upoly)
3821                         goto error;
3822         }
3823         qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3824         if (!qp->div)
3825                 goto error;
3826         for (i = 0; i < qp->div->n_row; ++i)
3827                 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3828
3829         isl_space_free(qp->dim);
3830         qp->dim = dim;
3831
3832         return qp;
3833 error:
3834         isl_space_free(dim);
3835         isl_qpolynomial_free(qp);
3836         return NULL;
3837 }
3838
3839 /* For each parameter or variable that does not appear in qp,
3840  * first eliminate the variable from all constraints and then set it to zero.
3841  */
3842 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3843         __isl_keep isl_qpolynomial *qp)
3844 {
3845         int *active = NULL;
3846         int i;
3847         int d;
3848         unsigned nparam;
3849         unsigned nvar;
3850
3851         if (!set || !qp)
3852                 goto error;
3853
3854         d = isl_space_dim(set->dim, isl_dim_all);
3855         active = isl_calloc_array(set->ctx, int, d);
3856         if (set_active(qp, active) < 0)
3857                 goto error;
3858
3859         for (i = 0; i < d; ++i)
3860                 if (!active[i])
3861                         break;
3862
3863         if (i == d) {
3864                 free(active);
3865                 return set;
3866         }
3867
3868         nparam = isl_space_dim(set->dim, isl_dim_param);
3869         nvar = isl_space_dim(set->dim, isl_dim_set);
3870         for (i = 0; i < nparam; ++i) {
3871                 if (active[i])
3872                         continue;
3873                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3874                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3875         }
3876         for (i = 0; i < nvar; ++i) {
3877                 if (active[nparam + i])
3878                         continue;
3879                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3880                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3881         }
3882
3883         free(active);
3884
3885         return set;
3886 error:
3887         free(active);
3888         isl_set_free(set);
3889         return NULL;
3890 }
3891
3892 struct isl_opt_data {
3893         isl_qpolynomial *qp;
3894         int first;
3895         isl_qpolynomial *opt;
3896         int max;
3897 };
3898
3899 static int opt_fn(__isl_take isl_point *pnt, void *user)
3900 {
3901         struct isl_opt_data *data = (struct isl_opt_data *)user;
3902         isl_qpolynomial *val;
3903
3904         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3905         if (data->first) {
3906                 data->first = 0;
3907                 data->opt = val;
3908         } else if (data->max) {
3909                 data->opt = isl_qpolynomial_max_cst(data->opt, val);
3910         } else {
3911                 data->opt = isl_qpolynomial_min_cst(data->opt, val);
3912         }
3913
3914         return 0;
3915 }
3916
3917 __isl_give isl_qpolynomial *isl_qpolynomial_opt_on_domain(
3918         __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3919 {
3920         struct isl_opt_data data = { NULL, 1, NULL, max };
3921
3922         if (!set || !qp)
3923                 goto error;
3924
3925         if (isl_upoly_is_cst(qp->upoly)) {
3926                 isl_set_free(set);
3927                 return qp;
3928         }
3929
3930         set = fix_inactive(set, qp);
3931
3932         data.qp = qp;
3933         if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3934                 goto error;
3935
3936         if (data.first) {
3937                 isl_space *space = isl_qpolynomial_get_domain_space(qp);
3938                 data.opt = isl_qpolynomial_zero_on_domain(space);
3939         }
3940
3941         isl_set_free(set);
3942         isl_qpolynomial_free(qp);
3943         return data.opt;
3944 error:
3945         isl_set_free(set);
3946         isl_qpolynomial_free(qp);
3947         isl_qpolynomial_free(data.opt);
3948         return NULL;
3949 }
3950
3951 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
3952         __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
3953 {
3954         int i;
3955         int n_sub;
3956         isl_ctx *ctx;
3957         struct isl_upoly **subs;
3958         isl_mat *mat, *diag;
3959
3960         qp = isl_qpolynomial_cow(qp);
3961         if (!qp || !morph)
3962                 goto error;
3963
3964         ctx = qp->dim->ctx;
3965         isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
3966
3967         n_sub = morph->inv->n_row - 1;
3968         if (morph->inv->n_row != morph->inv->n_col)
3969                 n_sub += qp->div->n_row;
3970         subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
3971         if (!subs)
3972                 goto error;
3973
3974         for (i = 0; 1 + i < morph->inv->n_row; ++i)
3975                 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
3976                                         morph->inv->row[0][0], morph->inv->n_col);
3977         if (morph->inv->n_row != morph->inv->n_col)
3978                 for (i = 0; i < qp->div->n_row; ++i)
3979                         subs[morph->inv->n_row - 1 + i] =
3980                             isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
3981
3982         qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
3983
3984         for (i = 0; i < n_sub; ++i)
3985                 isl_upoly_free(subs[i]);
3986         free(subs);
3987
3988         diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
3989         mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
3990         diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
3991         mat = isl_mat_diagonal(mat, diag);
3992         qp->div = isl_mat_product(qp->div, mat);
3993         isl_space_free(qp->dim);
3994         qp->dim = isl_space_copy(morph->ran->dim);
3995
3996         if (!qp->upoly || !qp->div || !qp->dim)
3997                 goto error;
3998
3999         isl_morph_free(morph);
4000
4001         return qp;
4002 error:
4003         isl_qpolynomial_free(qp);
4004         isl_morph_free(morph);
4005         return NULL;
4006 }
4007
4008 static int neg_entry(void **entry, void *user)
4009 {
4010         isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4011
4012         *pwqp = isl_pw_qpolynomial_neg(*pwqp);
4013
4014         return *pwqp ? 0 : -1;
4015 }
4016
4017 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
4018         __isl_take isl_union_pw_qpolynomial *upwqp)
4019 {
4020         upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4021         if (!upwqp)
4022                 return NULL;
4023
4024         if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4025                                    &neg_entry, NULL) < 0)
4026                 goto error;
4027
4028         return upwqp;
4029 error:
4030         isl_union_pw_qpolynomial_free(upwqp);
4031         return NULL;
4032 }
4033
4034 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4035         __isl_take isl_union_pw_qpolynomial *upwqp1,
4036         __isl_take isl_union_pw_qpolynomial *upwqp2)
4037 {
4038         return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
4039 }
4040
4041 /* Reorder the columns of the given div definitions according to the
4042  * given reordering.
4043  */
4044 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4045         __isl_take isl_reordering *r)
4046 {
4047         int i, j;
4048         isl_mat *mat;
4049         int extra;
4050
4051         if (!div || !r)
4052                 goto error;
4053
4054         extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4055         mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4056         if (!mat)
4057                 goto error;
4058
4059         for (i = 0; i < div->n_row; ++i) {
4060                 isl_seq_cpy(mat->row[i], div->row[i], 2);
4061                 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4062                 for (j = 0; j < r->len; ++j)
4063                         isl_int_set(mat->row[i][2 + r->pos[j]],
4064                                     div->row[i][2 + j]);
4065         }
4066
4067         isl_reordering_free(r);
4068         isl_mat_free(div);
4069         return mat;
4070 error:
4071         isl_reordering_free(r);
4072         isl_mat_free(div);
4073         return NULL;
4074 }
4075
4076 /* Reorder the dimension of "qp" according to the given reordering.
4077  */
4078 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4079         __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4080 {
4081         qp = isl_qpolynomial_cow(qp);
4082         if (!qp)
4083                 goto error;
4084
4085         r = isl_reordering_extend(r, qp->div->n_row);
4086         if (!r)
4087                 goto error;
4088
4089         qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4090         if (!qp->div)
4091                 goto error;
4092
4093         qp->upoly = reorder(qp->upoly, r->pos);
4094         if (!qp->upoly)
4095                 goto error;
4096
4097         qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4098
4099         isl_reordering_free(r);
4100         return qp;
4101 error:
4102         isl_qpolynomial_free(qp);
4103         isl_reordering_free(r);
4104         return NULL;
4105 }
4106
4107 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4108         __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4109 {
4110         if (!qp || !model)
4111                 goto error;
4112
4113         if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4114                 isl_reordering *exp;
4115
4116                 model = isl_space_drop_dims(model, isl_dim_in,
4117                                         0, isl_space_dim(model, isl_dim_in));
4118                 model = isl_space_drop_dims(model, isl_dim_out,
4119                                         0, isl_space_dim(model, isl_dim_out));
4120                 exp = isl_parameter_alignment_reordering(qp->dim, model);
4121                 exp = isl_reordering_extend_space(exp,
4122                                         isl_qpolynomial_get_domain_space(qp));
4123                 qp = isl_qpolynomial_realign_domain(qp, exp);
4124         }
4125
4126         isl_space_free(model);
4127         return qp;
4128 error:
4129         isl_space_free(model);
4130         isl_qpolynomial_free(qp);
4131         return NULL;
4132 }
4133
4134 struct isl_split_periods_data {
4135         int max_periods;
4136         isl_pw_qpolynomial *res;
4137 };
4138
4139 /* Create a slice where the integer division "div" has the fixed value "v".
4140  * In particular, if "div" refers to floor(f/m), then create a slice
4141  *
4142  *      m v <= f <= m v + (m - 1)
4143  *
4144  * or
4145  *
4146  *      f - m v >= 0
4147  *      -f + m v + (m - 1) >= 0
4148  */
4149 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4150         __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4151 {
4152         int total;
4153         isl_basic_set *bset = NULL;
4154         int k;
4155
4156         if (!dim || !qp)
4157                 goto error;
4158
4159         total = isl_space_dim(dim, isl_dim_all);
4160         bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4161
4162         k = isl_basic_set_alloc_inequality(bset);
4163         if (k < 0)
4164                 goto error;
4165         isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4166         isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4167
4168         k = isl_basic_set_alloc_inequality(bset);
4169         if (k < 0)
4170                 goto error;
4171         isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4172         isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4173         isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4174         isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4175
4176         isl_space_free(dim);
4177         return isl_set_from_basic_set(bset);
4178 error:
4179         isl_basic_set_free(bset);
4180         isl_space_free(dim);
4181         return NULL;
4182 }
4183
4184 static int split_periods(__isl_take isl_set *set,
4185         __isl_take isl_qpolynomial *qp, void *user);
4186
4187 /* Create a slice of the domain "set" such that integer division "div"
4188  * has the fixed value "v" and add the results to data->res,
4189  * replacing the integer division by "v" in "qp".
4190  */
4191 static int set_div(__isl_take isl_set *set,
4192         __isl_take isl_qpolynomial *qp, int div, isl_int v,
4193         struct isl_split_periods_data *data)
4194 {
4195         int i;
4196         int total;
4197         isl_set *slice;
4198         struct isl_upoly *cst;
4199
4200         slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4201         set = isl_set_intersect(set, slice);
4202
4203         if (!qp)
4204                 goto error;
4205
4206         total = isl_space_dim(qp->dim, isl_dim_all);
4207
4208         for (i = div + 1; i < qp->div->n_row; ++i) {
4209                 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4210                         continue;
4211                 isl_int_addmul(qp->div->row[i][1],
4212                                 qp->div->row[i][2 + total + div], v);
4213                 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4214         }
4215
4216         cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4217         qp = substitute_div(qp, div, cst);
4218
4219         return split_periods(set, qp, data);
4220 error:
4221         isl_set_free(set);
4222         isl_qpolynomial_free(qp);
4223         return -1;
4224 }
4225
4226 /* Split the domain "set" such that integer division "div"
4227  * has a fixed value (ranging from "min" to "max") on each slice
4228  * and add the results to data->res.
4229  */
4230 static int split_div(__isl_take isl_set *set,
4231         __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4232         struct isl_split_periods_data *data)
4233 {
4234         for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4235                 isl_set *set_i = isl_set_copy(set);
4236                 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4237
4238                 if (set_div(set_i, qp_i, div, min, data) < 0)
4239                         goto error;
4240         }
4241         isl_set_free(set);
4242         isl_qpolynomial_free(qp);
4243         return 0;
4244 error:
4245         isl_set_free(set);
4246         isl_qpolynomial_free(qp);
4247         return -1;
4248 }
4249
4250 /* If "qp" refers to any integer division
4251  * that can only attain "max_periods" distinct values on "set"
4252  * then split the domain along those distinct values.
4253  * Add the results (or the original if no splitting occurs)
4254  * to data->res.
4255  */
4256 static int split_periods(__isl_take isl_set *set,
4257         __isl_take isl_qpolynomial *qp, void *user)
4258 {
4259         int i;
4260         isl_pw_qpolynomial *pwqp;
4261         struct isl_split_periods_data *data;
4262         isl_int min, max;
4263         int total;
4264         int r = 0;
4265
4266         data = (struct isl_split_periods_data *)user;
4267
4268         if (!set || !qp)
4269                 goto error;
4270
4271         if (qp->div->n_row == 0) {
4272                 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4273                 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4274                 return 0;
4275         }
4276
4277         isl_int_init(min);
4278         isl_int_init(max);
4279         total = isl_space_dim(qp->dim, isl_dim_all);
4280         for (i = 0; i < qp->div->n_row; ++i) {
4281                 enum isl_lp_result lp_res;
4282
4283                 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4284                                                 qp->div->n_row) != -1)
4285                         continue;
4286
4287                 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4288                                           set->ctx->one, &min, NULL, NULL);
4289                 if (lp_res == isl_lp_error)
4290                         goto error2;
4291                 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4292                         continue;
4293                 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4294
4295                 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4296                                           set->ctx->one, &max, NULL, NULL);
4297                 if (lp_res == isl_lp_error)
4298                         goto error2;
4299                 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4300                         continue;
4301                 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4302
4303                 isl_int_sub(max, max, min);
4304                 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4305                         isl_int_add(max, max, min);
4306                         break;
4307                 }
4308         }
4309
4310         if (i < qp->div->n_row) {
4311                 r = split_div(set, qp, i, min, max, data);
4312         } else {
4313                 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4314                 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4315         }
4316
4317         isl_int_clear(max);
4318         isl_int_clear(min);
4319
4320         return r;
4321 error2:
4322         isl_int_clear(max);
4323         isl_int_clear(min);
4324 error:
4325         isl_set_free(set);
4326         isl_qpolynomial_free(qp);
4327         return -1;
4328 }
4329
4330 /* If any quasi-polynomial in pwqp refers to any integer division
4331  * that can only attain "max_periods" distinct values on its domain
4332  * then split the domain along those distinct values.
4333  */
4334 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4335         __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4336 {
4337         struct isl_split_periods_data data;
4338
4339         data.max_periods = max_periods;
4340         data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4341
4342         if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4343                 goto error;
4344
4345         isl_pw_qpolynomial_free(pwqp);
4346
4347         return data.res;
4348 error:
4349         isl_pw_qpolynomial_free(data.res);
4350         isl_pw_qpolynomial_free(pwqp);
4351         return NULL;
4352 }
4353
4354 /* Construct a piecewise quasipolynomial that is constant on the given
4355  * domain.  In particular, it is
4356  *      0       if cst == 0
4357  *      1       if cst == 1
4358  *  infinity    if cst == -1
4359  */
4360 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4361         __isl_take isl_basic_set *bset, int cst)
4362 {
4363         isl_space *dim;
4364         isl_qpolynomial *qp;
4365
4366         if (!bset)
4367                 return NULL;
4368
4369         bset = isl_basic_set_params(bset);
4370         dim = isl_basic_set_get_space(bset);
4371         if (cst < 0)
4372                 qp = isl_qpolynomial_infty_on_domain(dim);
4373         else if (cst == 0)
4374                 qp = isl_qpolynomial_zero_on_domain(dim);
4375         else
4376                 qp = isl_qpolynomial_one_on_domain(dim);
4377         return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4378 }
4379
4380 /* Factor bset, call fn on each of the factors and return the product.
4381  *
4382  * If no factors can be found, simply call fn on the input.
4383  * Otherwise, construct the factors based on the factorizer,
4384  * call fn on each factor and compute the product.
4385  */
4386 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4387         __isl_take isl_basic_set *bset,
4388         __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4389 {
4390         int i, n;
4391         isl_space *dim;
4392         isl_set *set;
4393         isl_factorizer *f;
4394         isl_qpolynomial *qp;
4395         isl_pw_qpolynomial *pwqp;
4396         unsigned nparam;
4397         unsigned nvar;
4398
4399         f = isl_basic_set_factorizer(bset);
4400         if (!f)
4401                 goto error;
4402         if (f->n_group == 0) {
4403                 isl_factorizer_free(f);
4404                 return fn(bset);
4405         }
4406
4407         nparam = isl_basic_set_dim(bset, isl_dim_param);
4408         nvar = isl_basic_set_dim(bset, isl_dim_set);
4409
4410         dim = isl_basic_set_get_space(bset);
4411         dim = isl_space_domain(dim);
4412         set = isl_set_universe(isl_space_copy(dim));
4413         qp = isl_qpolynomial_one_on_domain(dim);
4414         pwqp = isl_pw_qpolynomial_alloc(set, qp);
4415
4416         bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4417
4418         for (i = 0, n = 0; i < f->n_group; ++i) {
4419                 isl_basic_set *bset_i;
4420                 isl_pw_qpolynomial *pwqp_i;
4421
4422                 bset_i = isl_basic_set_copy(bset);
4423                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4424                             nparam + n + f->len[i], nvar - n - f->len[i]);
4425                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4426                             nparam, n);
4427                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4428                             n + f->len[i], nvar - n - f->len[i]);
4429                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4430
4431                 pwqp_i = fn(bset_i);
4432                 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4433
4434                 n += f->len[i];
4435         }
4436
4437         isl_basic_set_free(bset);
4438         isl_factorizer_free(f);
4439
4440         return pwqp;
4441 error:
4442         isl_basic_set_free(bset);
4443         return NULL;
4444 }
4445
4446 /* Factor bset, call fn on each of the factors and return the product.
4447  * The function is assumed to evaluate to zero on empty domains,
4448  * to one on zero-dimensional domains and to infinity on unbounded domains
4449  * and will not be called explicitly on zero-dimensional or unbounded domains.
4450  *
4451  * We first check for some special cases and remove all equalities.
4452  * Then we hand over control to compressed_multiplicative_call.
4453  */
4454 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4455         __isl_take isl_basic_set *bset,
4456         __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4457 {
4458         int bounded;
4459         isl_morph *morph;
4460         isl_pw_qpolynomial *pwqp;
4461
4462         if (!bset)
4463                 return NULL;
4464
4465         if (isl_basic_set_plain_is_empty(bset))
4466                 return constant_on_domain(bset, 0);
4467
4468         if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4469                 return constant_on_domain(bset, 1);
4470
4471         bounded = isl_basic_set_is_bounded(bset);
4472         if (bounded < 0)
4473                 goto error;
4474         if (!bounded)
4475                 return constant_on_domain(bset, -1);
4476
4477         if (bset->n_eq == 0)
4478                 return compressed_multiplicative_call(bset, fn);
4479
4480         morph = isl_basic_set_full_compression(bset);
4481         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4482
4483         pwqp = compressed_multiplicative_call(bset, fn);
4484
4485         morph = isl_morph_dom_params(morph);
4486         morph = isl_morph_ran_params(morph);
4487         morph = isl_morph_inverse(morph);
4488
4489         pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4490
4491         return pwqp;
4492 error:
4493         isl_basic_set_free(bset);
4494         return NULL;
4495 }
4496
4497 /* Drop all floors in "qp", turning each integer division [a/m] into
4498  * a rational division a/m.  If "down" is set, then the integer division
4499  * is replaces by (a-(m-1))/m instead.
4500  */
4501 static __isl_give isl_qpolynomial *qp_drop_floors(
4502         __isl_take isl_qpolynomial *qp, int down)
4503 {
4504         int i;
4505         struct isl_upoly *s;
4506
4507         if (!qp)
4508                 return NULL;
4509         if (qp->div->n_row == 0)
4510                 return qp;
4511
4512         qp = isl_qpolynomial_cow(qp);
4513         if (!qp)
4514                 return NULL;
4515
4516         for (i = qp->div->n_row - 1; i >= 0; --i) {
4517                 if (down) {
4518                         isl_int_sub(qp->div->row[i][1],
4519                                     qp->div->row[i][1], qp->div->row[i][0]);
4520                         isl_int_add_ui(qp->div->row[i][1],
4521                                        qp->div->row[i][1], 1);
4522                 }
4523                 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4524                                         qp->div->row[i][0], qp->div->n_col - 1);
4525                 qp = substitute_div(qp, i, s);
4526                 if (!qp)
4527                         return NULL;
4528         }
4529
4530         return qp;
4531 }
4532
4533 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4534  * a rational division a/m.
4535  */
4536 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4537         __isl_take isl_pw_qpolynomial *pwqp)
4538 {
4539         int i;
4540
4541         if (!pwqp)
4542                 return NULL;
4543
4544         if (isl_pw_qpolynomial_is_zero(pwqp))
4545                 return pwqp;
4546
4547         pwqp = isl_pw_qpolynomial_cow(pwqp);
4548         if (!pwqp)
4549                 return NULL;
4550
4551         for (i = 0; i < pwqp->n; ++i) {
4552                 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4553                 if (!pwqp->p[i].qp)
4554                         goto error;
4555         }
4556
4557         return pwqp;
4558 error:
4559         isl_pw_qpolynomial_free(pwqp);
4560         return NULL;
4561 }
4562
4563 /* Adjust all the integer divisions in "qp" such that they are at least
4564  * one over the given orthant (identified by "signs").  This ensures
4565  * that they will still be non-negative even after subtracting (m-1)/m.
4566  *
4567  * In particular, f is replaced by f' + v, changing f = [a/m]
4568  * to f' = [(a - m v)/m].
4569  * If the constant term k in a is smaller than m,
4570  * the constant term of v is set to floor(k/m) - 1.
4571  * For any other term, if the coefficient c and the variable x have
4572  * the same sign, then no changes are needed.
4573  * Otherwise, if the variable is positive (and c is negative),
4574  * then the coefficient of x in v is set to floor(c/m).
4575  * If the variable is negative (and c is positive),
4576  * then the coefficient of x in v is set to ceil(c/m).
4577  */
4578 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4579         int *signs)
4580 {
4581         int i, j;
4582         int total;
4583         isl_vec *v = NULL;
4584         struct isl_upoly *s;
4585
4586         qp = isl_qpolynomial_cow(qp);
4587         if (!qp)
4588                 return NULL;
4589         qp->div = isl_mat_cow(qp->div);
4590         if (!qp->div)
4591                 goto error;
4592
4593         total = isl_space_dim(qp->dim, isl_dim_all);
4594         v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4595
4596         for (i = 0; i < qp->div->n_row; ++i) {
4597                 isl_int *row = qp->div->row[i];
4598                 v = isl_vec_clr(v);
4599                 if (!v)
4600                         goto error;
4601                 if (isl_int_lt(row[1], row[0])) {
4602                         isl_int_fdiv_q(v->el[0], row[1], row[0]);
4603                         isl_int_sub_ui(v->el[0], v->el[0], 1);
4604                         isl_int_submul(row[1], row[0], v->el[0]);
4605                 }
4606                 for (j = 0; j < total; ++j) {
4607                         if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4608                                 continue;
4609                         if (signs[j] < 0)
4610                                 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4611                         else
4612                                 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4613                         isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4614                 }
4615                 for (j = 0; j < i; ++j) {
4616                         if (isl_int_sgn(row[2 + total + j]) >= 0)
4617                                 continue;
4618                         isl_int_fdiv_q(v->el[1 + total + j],
4619                                         row[2 + total + j], row[0]);
4620                         isl_int_submul(row[2 + total + j],
4621                                         row[0], v->el[1 + total + j]);
4622                 }
4623                 for (j = i + 1; j < qp->div->n_row; ++j) {
4624                         if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4625                                 continue;
4626                         isl_seq_combine(qp->div->row[j] + 1,
4627                                 qp->div->ctx->one, qp->div->row[j] + 1,
4628                                 qp->div->row[j][2 + total + i], v->el, v->size);
4629                 }
4630                 isl_int_set_si(v->el[1 + total + i], 1);
4631                 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4632                                         qp->div->ctx->one, v->size);
4633                 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4634                 isl_upoly_free(s);
4635                 if (!qp->upoly)
4636                         goto error;
4637         }
4638
4639         isl_vec_free(v);
4640         return qp;
4641 error:
4642         isl_vec_free(v);
4643         isl_qpolynomial_free(qp);
4644         return NULL;
4645 }
4646
4647 struct isl_to_poly_data {
4648         int sign;
4649         isl_pw_qpolynomial *res;
4650         isl_qpolynomial *qp;
4651 };
4652
4653 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4654  * We first make all integer divisions positive and then split the
4655  * quasipolynomials into terms with sign data->sign (the direction
4656  * of the requested approximation) and terms with the opposite sign.
4657  * In the first set of terms, each integer division [a/m] is
4658  * overapproximated by a/m, while in the second it is underapproximated
4659  * by (a-(m-1))/m.
4660  */
4661 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4662         void *user)
4663 {
4664         struct isl_to_poly_data *data = user;
4665         isl_pw_qpolynomial *t;
4666         isl_qpolynomial *qp, *up, *down;
4667
4668         qp = isl_qpolynomial_copy(data->qp);
4669         qp = make_divs_pos(qp, signs);
4670
4671         up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4672         up = qp_drop_floors(up, 0);
4673         down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4674         down = qp_drop_floors(down, 1);
4675
4676         isl_qpolynomial_free(qp);
4677         qp = isl_qpolynomial_add(up, down);
4678
4679         t = isl_pw_qpolynomial_alloc(orthant, qp);
4680         data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4681
4682         return 0;
4683 }
4684
4685 /* Approximate each quasipolynomial by a polynomial.  If "sign" is positive,
4686  * the polynomial will be an overapproximation.  If "sign" is negative,
4687  * it will be an underapproximation.  If "sign" is zero, the approximation
4688  * will lie somewhere in between.
4689  *
4690  * In particular, is sign == 0, we simply drop the floors, turning
4691  * the integer divisions into rational divisions.
4692  * Otherwise, we split the domains into orthants, make all integer divisions
4693  * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4694  * depending on the requested sign and the sign of the term in which
4695  * the integer division appears.
4696  */
4697 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4698         __isl_take isl_pw_qpolynomial *pwqp, int sign)
4699 {
4700         int i;
4701         struct isl_to_poly_data data;
4702
4703         if (sign == 0)
4704                 return pwqp_drop_floors(pwqp);
4705
4706         if (!pwqp)
4707                 return NULL;
4708
4709         data.sign = sign;
4710         data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4711
4712         for (i = 0; i < pwqp->n; ++i) {
4713                 if (pwqp->p[i].qp->div->n_row == 0) {
4714                         isl_pw_qpolynomial *t;
4715                         t = isl_pw_qpolynomial_alloc(
4716                                         isl_set_copy(pwqp->p[i].set),
4717                                         isl_qpolynomial_copy(pwqp->p[i].qp));
4718                         data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4719                         continue;
4720                 }
4721                 data.qp = pwqp->p[i].qp;
4722                 if (isl_set_foreach_orthant(pwqp->p[i].set,
4723                                         &to_polynomial_on_orthant, &data) < 0)
4724                         goto error;
4725         }
4726
4727         isl_pw_qpolynomial_free(pwqp);
4728
4729         return data.res;
4730 error:
4731         isl_pw_qpolynomial_free(pwqp);
4732         isl_pw_qpolynomial_free(data.res);
4733         return NULL;
4734 }
4735
4736 static int poly_entry(void **entry, void *user)
4737 {
4738         int *sign = user;
4739         isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4740
4741         *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4742
4743         return *pwqp ? 0 : -1;
4744 }
4745
4746 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4747         __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4748 {
4749         upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4750         if (!upwqp)
4751                 return NULL;
4752
4753         if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4754                                    &poly_entry, &sign) < 0)
4755                 goto error;
4756
4757         return upwqp;
4758 error:
4759         isl_union_pw_qpolynomial_free(upwqp);
4760         return NULL;
4761 }
4762
4763 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4764         __isl_take isl_qpolynomial *qp)
4765 {
4766         int i, k;
4767         isl_space *dim;
4768         isl_vec *aff = NULL;
4769         isl_basic_map *bmap = NULL;
4770         unsigned pos;
4771         unsigned n_div;
4772
4773         if (!qp)
4774                 return NULL;
4775         if (!isl_upoly_is_affine(qp->upoly))
4776                 isl_die(qp->dim->ctx, isl_error_invalid,
4777                         "input quasi-polynomial not affine", goto error);
4778         aff = isl_qpolynomial_extract_affine(qp);
4779         if (!aff)
4780                 goto error;
4781         dim = isl_qpolynomial_get_space(qp);
4782         pos = 1 + isl_space_offset(dim, isl_dim_out);
4783         n_div = qp->div->n_row;
4784         bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4785
4786         for (i = 0; i < n_div; ++i) {
4787                 k = isl_basic_map_alloc_div(bmap);
4788                 if (k < 0)
4789                         goto error;
4790                 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4791                 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4792                 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4793                         goto error;
4794         }
4795         k = isl_basic_map_alloc_equality(bmap);
4796         if (k < 0)
4797                 goto error;
4798         isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4799         isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4800         isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4801
4802         isl_vec_free(aff);
4803         isl_qpolynomial_free(qp);
4804         bmap = isl_basic_map_finalize(bmap);
4805         return bmap;
4806 error:
4807         isl_vec_free(aff);
4808         isl_qpolynomial_free(qp);
4809         isl_basic_map_free(bmap);
4810         return NULL;
4811 }