isl_polynomial.c

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