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