isl_polynomial.c

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