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