isl_polynomial.c

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