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