isl_polynomial.c

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