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