isl_polynomial.c

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