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