isl_polynomial.c

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