isl_polynomial.c

   1 #include <stdlib.h>
   2 #include <isl_seq.h>
   3 #include <isl_polynomial_private.h>
   4 #include <isl_point_private.h>
   5 #include <isl_dim_private.h>
   6 #include <isl_map_private.h>
   7 #include <isl_map_private.h>
   8
   9 static unsigned pos(__isl_keep isl_dim *dim, enum isl_dim_type type)
  10 {
  11         switch (type) {
  12         case isl_dim_param:     return 0;
  13         case isl_dim_in:        return dim->nparam;
  14         case isl_dim_out:       return dim->nparam + dim->n_in;
  15         }
  16 }
  17
  18 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
  19 {
  20         if (!up)
  21                 return -1;
  22
  23         return up->var < 0;
  24 }
  25
  26 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
  27 {
  28         if (!up)
  29                 return NULL;
  30
  31         isl_assert(up->ctx, up->var < 0, return NULL);
  32
  33         return (struct isl_upoly_cst *)up;
  34 }
  35
  36 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
  37 {
  38         if (!up)
  39                 return NULL;
  40
  41         isl_assert(up->ctx, up->var >= 0, return NULL);
  42
  43         return (struct isl_upoly_rec *)up;
  44 }
  45
  46 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
  47         __isl_keep struct isl_upoly *up2)
  48 {
  49         int i;
  50         struct isl_upoly_rec *rec1, *rec2;
  51
  52         if (!up1 || !up2)
  53                 return -1;
  54         if (up1 == up2)
  55                 return 1;
  56         if (up1->var != up2->var)
  57                 return 0;
  58         if (isl_upoly_is_cst(up1)) {
  59                 struct isl_upoly_cst *cst1, *cst2;
  60                 cst1 = isl_upoly_as_cst(up1);
  61                 cst2 = isl_upoly_as_cst(up2);
  62                 if (!cst1 || !cst2)
  63                         return -1;
  64                 return isl_int_eq(cst1->n, cst2->n) &&
  65                        isl_int_eq(cst1->d, cst2->d);
  66         }
  67
  68         rec1 = isl_upoly_as_rec(up1);
  69         rec2 = isl_upoly_as_rec(up2);
  70         if (!rec1 || !rec2)
  71                 return -1;
  72
  73         if (rec1->n != rec2->n)
  74                 return 0;
  75
  76         for (i = 0; i < rec1->n; ++i) {
  77                 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
  78                 if (eq < 0 || !eq)
  79                         return eq;
  80         }
  81
  82         return 1;
  83 }
  84
  85 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
  86 {
  87         struct isl_upoly_cst *cst;
  88
  89         if (!up)
  90                 return -1;
  91         if (!isl_upoly_is_cst(up))
  92                 return 0;
  93
  94         cst = isl_upoly_as_cst(up);
  95         if (!cst)
  96                 return -1;
  97
  98         return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
  99 }
 100
 101 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
 102 {
 103         struct isl_upoly_cst *cst;
 104
 105         if (!up)
 106                 return -1;
 107         if (!isl_upoly_is_cst(up))
 108                 return 0;
 109
 110         cst = isl_upoly_as_cst(up);
 111         if (!cst)
 112                 return -1;
 113
 114         return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
 115 }
 116
 117 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
 118 {
 119         struct isl_upoly_cst *cst;
 120
 121         if (!up)
 122                 return -1;
 123         if (!isl_upoly_is_cst(up))
 124                 return 0;
 125
 126         cst = isl_upoly_as_cst(up);
 127         if (!cst)
 128                 return -1;
 129
 130         return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
 131 }
 132
 133 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
 134 {
 135         struct isl_upoly_cst *cst;
 136
 137         if (!up)
 138                 return -1;
 139         if (!isl_upoly_is_cst(up))
 140                 return 0;
 141
 142         cst = isl_upoly_as_cst(up);
 143         if (!cst)
 144                 return -1;
 145
 146         return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
 147 }
 148
 149 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
 150 {
 151         struct isl_upoly_cst *cst;
 152
 153         if (!up)
 154                 return -1;
 155         if (!isl_upoly_is_cst(up))
 156                 return 0;
 157
 158         cst = isl_upoly_as_cst(up);
 159         if (!cst)
 160                 return -1;
 161
 162         return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
 163 }
 164
 165 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
 166 {
 167         struct isl_upoly_cst *cst;
 168
 169         if (!up)
 170                 return -1;
 171         if (!isl_upoly_is_cst(up))
 172                 return 0;
 173
 174         cst = isl_upoly_as_cst(up);
 175         if (!cst)
 176                 return -1;
 177
 178         return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
 179 }
 180
 181 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
 182 {
 183         struct isl_upoly_cst *cst;
 184
 185         cst = isl_alloc_type(ctx, struct isl_upoly_cst);
 186         if (!cst)
 187                 return NULL;
 188
 189         cst->up.ref = 1;
 190         cst->up.ctx = ctx;
 191         isl_ctx_ref(ctx);
 192         cst->up.var = -1;
 193
 194         isl_int_init(cst->n);
 195         isl_int_init(cst->d);
 196
 197         return cst;
 198 }
 199
 200 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
 201 {
 202         struct isl_upoly_cst *cst;
 203
 204         cst = isl_upoly_cst_alloc(ctx);
 205         if (!cst)
 206                 return NULL;
 207
 208         isl_int_set_si(cst->n, 0);
 209         isl_int_set_si(cst->d, 1);
 210
 211         return &cst->up;
 212 }
 213
 214 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
 215 {
 216         struct isl_upoly_cst *cst;
 217
 218         cst = isl_upoly_cst_alloc(ctx);
 219         if (!cst)
 220                 return NULL;
 221
 222         isl_int_set_si(cst->n, 1);
 223         isl_int_set_si(cst->d, 0);
 224
 225         return &cst->up;
 226 }
 227
 228 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
 229 {
 230         struct isl_upoly_cst *cst;
 231
 232         cst = isl_upoly_cst_alloc(ctx);
 233         if (!cst)
 234                 return NULL;
 235
 236         isl_int_set_si(cst->n, 0);
 237         isl_int_set_si(cst->d, 0);
 238
 239         return &cst->up;
 240 }
 241
 242 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
 243         isl_int n, isl_int d)
 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(cst->n, n);
 252         isl_int_set(cst->d, d);
 253
 254         return &cst->up;
 255 }
 256
 257 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
 258         int var, int size)
 259 {
 260         struct isl_upoly_rec *rec;
 261
 262         isl_assert(ctx, var >= 0, return NULL);
 263         isl_assert(ctx, size >= 0, return NULL);
 264         rec = isl_calloc(dim->ctx, struct isl_upoly_rec,
 265                         sizeof(struct isl_upoly_rec) +
 266                         (size - 1) * sizeof(struct isl_upoly *));
 267         if (!rec)
 268                 return NULL;
 269
 270         rec->up.ref = 1;
 271         rec->up.ctx = ctx;
 272         isl_ctx_ref(ctx);
 273         rec->up.var = var;
 274
 275         rec->n = 0;
 276         rec->size = size;
 277
 278         return rec;
 279 error:
 280         isl_upoly_free(&rec->up);
 281         return NULL;
 282 }
 283
 284 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
 285 {
 286         struct isl_upoly_cst *cst;
 287
 288         if (!qp)
 289                 return -1;
 290
 291         return isl_upoly_is_zero(qp->upoly);
 292 }
 293
 294 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
 295 {
 296         struct isl_upoly_cst *cst;
 297
 298         if (!qp)
 299                 return -1;
 300
 301         return isl_upoly_is_one(qp->upoly);
 302 }
 303
 304 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
 305 {
 306         isl_int_clear(cst->n);
 307         isl_int_clear(cst->d);
 308 }
 309
 310 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
 311 {
 312         int i;
 313
 314         for (i = 0; i < rec->n; ++i)
 315                 isl_upoly_free(rec->p[i]);
 316 }
 317
 318 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
 319 {
 320         if (!up)
 321                 return NULL;
 322
 323         up->ref++;
 324         return up;
 325 }
 326
 327 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
 328 {
 329         struct isl_upoly_cst *cst;
 330         struct isl_upoly_cst *dup;
 331
 332         cst = isl_upoly_as_cst(up);
 333         if (!cst)
 334                 return NULL;
 335
 336         dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
 337         if (!dup)
 338                 return NULL;
 339         isl_int_set(dup->n, cst->n);
 340         isl_int_set(dup->d, cst->d);
 341
 342         return &dup->up;
 343 }
 344
 345 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
 346 {
 347         int i;
 348         struct isl_upoly_rec *rec;
 349         struct isl_upoly_rec *dup;
 350
 351         rec = isl_upoly_as_rec(up);
 352         if (!rec)
 353                 return NULL;
 354
 355         dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
 356         if (!dup)
 357                 return NULL;
 358
 359         for (i = 0; i < rec->n; ++i) {
 360                 dup->p[i] = isl_upoly_copy(rec->p[i]);
 361                 if (!dup->p[i])
 362                         goto error;
 363                 dup->n++;
 364         }
 365
 366         return &dup->up;
 367 error:
 368         isl_upoly_free(&dup->up);
 369         return NULL;
 370 }
 371
 372 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
 373 {
 374         struct isl_upoly *dup;
 375
 376         if (!up)
 377                 return NULL;
 378
 379         if (isl_upoly_is_cst(up))
 380                 return isl_upoly_dup_cst(up);
 381         else
 382                 return isl_upoly_dup_rec(up);
 383 }
 384
 385 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
 386 {
 387         if (!up)
 388                 return NULL;
 389
 390         if (up->ref == 1)
 391                 return up;
 392         up->ref--;
 393         return isl_upoly_dup(up);
 394 }
 395
 396 void isl_upoly_free(__isl_take struct isl_upoly *up)
 397 {
 398         if (!up)
 399                 return;
 400
 401         if (--up->ref > 0)
 402                 return;
 403
 404         if (up->var < 0)
 405                 upoly_free_cst((struct isl_upoly_cst *)up);
 406         else
 407                 upoly_free_rec((struct isl_upoly_rec *)up);
 408
 409         isl_ctx_deref(up->ctx);
 410         free(up);
 411 }
 412
 413 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
 414 {
 415         isl_int gcd;
 416
 417         isl_int_init(gcd);
 418         isl_int_gcd(gcd, cst->n, cst->d);
 419         if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
 420                 isl_int_divexact(cst->n, cst->n, gcd);
 421                 isl_int_divexact(cst->d, cst->d, gcd);
 422         }
 423         isl_int_clear(gcd);
 424 }
 425
 426 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
 427         __isl_take struct isl_upoly *up2)
 428 {
 429         struct isl_upoly_cst *cst1;
 430         struct isl_upoly_cst *cst2;
 431
 432         up1 = isl_upoly_cow(up1);
 433         if (!up1 || !up2)
 434                 goto error;
 435
 436         cst1 = isl_upoly_as_cst(up1);
 437         cst2 = isl_upoly_as_cst(up2);
 438
 439         if (isl_int_eq(cst1->d, cst2->d))
 440                 isl_int_add(cst1->n, cst1->n, cst2->n);
 441         else {
 442                 isl_int_mul(cst1->n, cst1->n, cst2->d);
 443                 isl_int_addmul(cst1->n, cst2->n, cst1->d);
 444                 isl_int_mul(cst1->d, cst1->d, cst2->d);
 445         }
 446
 447         isl_upoly_cst_reduce(cst1);
 448
 449         isl_upoly_free(up2);
 450         return up1;
 451 error:
 452         isl_upoly_free(up1);
 453         isl_upoly_free(up2);
 454         return NULL;
 455 }
 456
 457 static __isl_give struct isl_upoly *replace_by_zero(
 458         __isl_take struct isl_upoly *up)
 459 {
 460         struct isl_ctx *ctx;
 461
 462         if (!up)
 463                 return NULL;
 464         ctx = up->ctx;
 465         isl_upoly_free(up);
 466         return isl_upoly_zero(ctx);
 467 }
 468
 469 static __isl_give struct isl_upoly *replace_by_constant_term(
 470         __isl_take struct isl_upoly *up)
 471 {
 472         struct isl_upoly_rec *rec;
 473         struct isl_upoly *cst;
 474
 475         if (!up)
 476                 return NULL;
 477
 478         rec = isl_upoly_as_rec(up);
 479         if (!rec)
 480                 goto error;
 481         cst = isl_upoly_copy(rec->p[0]);
 482         isl_upoly_free(up);
 483         return cst;
 484 error:
 485         isl_upoly_free(up);
 486         return NULL;
 487 }
 488
 489 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
 490         __isl_take struct isl_upoly *up2)
 491 {
 492         int i;
 493         struct isl_upoly_rec *rec1, *rec2;
 494
 495         if (!up1 || !up2)
 496                 goto error;
 497
 498         if (isl_upoly_is_nan(up1)) {
 499                 isl_upoly_free(up2);
 500                 return up1;
 501         }
 502
 503         if (isl_upoly_is_nan(up2)) {
 504                 isl_upoly_free(up1);
 505                 return up2;
 506         }
 507
 508         if (isl_upoly_is_zero(up1)) {
 509                 isl_upoly_free(up1);
 510                 return up2;
 511         }
 512
 513         if (isl_upoly_is_zero(up2)) {
 514                 isl_upoly_free(up2);
 515                 return up1;
 516         }
 517
 518         if (up1->var < up2->var)
 519                 return isl_upoly_sum(up2, up1);
 520
 521         if (up2->var < up1->var) {
 522                 struct isl_upoly_rec *rec;
 523                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 524                         isl_upoly_free(up1);
 525                         return up2;
 526                 }
 527                 up1 = isl_upoly_cow(up1);
 528                 rec = isl_upoly_as_rec(up1);
 529                 if (!rec)
 530                         goto error;
 531                 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
 532                 if (rec->n == 1)
 533                         up1 = replace_by_constant_term(up1);
 534                 return up1;
 535         }
 536
 537         if (isl_upoly_is_cst(up1))
 538                 return isl_upoly_sum_cst(up1, up2);
 539
 540         rec1 = isl_upoly_as_rec(up1);
 541         rec2 = isl_upoly_as_rec(up2);
 542         if (!rec1 || !rec2)
 543                 goto error;
 544
 545         if (rec1->n < rec2->n)
 546                 return isl_upoly_sum(up2, up1);
 547
 548         up1 = isl_upoly_cow(up1);
 549         rec1 = isl_upoly_as_rec(up1);
 550         if (!rec1)
 551                 goto error;
 552
 553         for (i = rec2->n - 1; i >= 0; --i) {
 554                 rec1->p[i] = isl_upoly_sum(rec1->p[i],
 555                                             isl_upoly_copy(rec2->p[i]));
 556                 if (!rec1->p[i])
 557                         goto error;
 558                 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
 559                         isl_upoly_free(rec1->p[i]);
 560                         rec1->n--;
 561                 }
 562         }
 563
 564         if (rec1->n == 0)
 565                 up1 = replace_by_zero(up1);
 566         else if (rec1->n == 1)
 567                 up1 = replace_by_constant_term(up1);
 568
 569         isl_upoly_free(up2);
 570
 571         return up1;
 572 error:
 573         isl_upoly_free(up1);
 574         isl_upoly_free(up2);
 575         return NULL;
 576 }
 577
 578 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
 579 {
 580         struct isl_upoly_cst *cst;
 581
 582         if (isl_upoly_is_zero(up))
 583                 return up;
 584
 585         up = isl_upoly_cow(up);
 586         if (!up)
 587                 return NULL;
 588
 589         cst = isl_upoly_as_cst(up);
 590
 591         isl_int_neg(cst->n, cst->n);
 592
 593         return up;
 594 }
 595
 596 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
 597 {
 598         int i;
 599         struct isl_upoly_rec *rec;
 600
 601         if (!up)
 602                 return NULL;
 603
 604         if (isl_upoly_is_cst(up))
 605                 return isl_upoly_neg_cst(up);
 606
 607         up = isl_upoly_cow(up);
 608         rec = isl_upoly_as_rec(up);
 609         if (!rec)
 610                 goto error;
 611
 612         for (i = 0; i < rec->n; ++i) {
 613                 rec->p[i] = isl_upoly_neg(rec->p[i]);
 614                 if (!rec->p[i])
 615                         goto error;
 616         }
 617
 618         return up;
 619 error:
 620         isl_upoly_free(up);
 621         return NULL;
 622 }
 623
 624 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
 625         __isl_take struct isl_upoly *up2)
 626 {
 627         struct isl_upoly_cst *cst1;
 628         struct isl_upoly_cst *cst2;
 629
 630         up1 = isl_upoly_cow(up1);
 631         if (!up1 || !up2)
 632                 goto error;
 633
 634         cst1 = isl_upoly_as_cst(up1);
 635         cst2 = isl_upoly_as_cst(up2);
 636
 637         isl_int_mul(cst1->n, cst1->n, cst2->n);
 638         isl_int_mul(cst1->d, cst1->d, cst2->d);
 639
 640         isl_upoly_cst_reduce(cst1);
 641
 642         isl_upoly_free(up2);
 643         return up1;
 644 error:
 645         isl_upoly_free(up1);
 646         isl_upoly_free(up2);
 647         return NULL;
 648 }
 649
 650 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
 651         __isl_take struct isl_upoly *up2)
 652 {
 653         struct isl_upoly_rec *rec1;
 654         struct isl_upoly_rec *rec2;
 655         struct isl_upoly_rec *res;
 656         int i, j;
 657         int size;
 658
 659         rec1 = isl_upoly_as_rec(up1);
 660         rec2 = isl_upoly_as_rec(up2);
 661         if (!rec1 || !rec2)
 662                 goto error;
 663         size = rec1->n + rec2->n - 1;
 664         res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
 665         if (!res)
 666                 goto error;
 667
 668         for (i = 0; i < rec1->n; ++i) {
 669                 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
 670                                             isl_upoly_copy(rec1->p[i]));
 671                 if (!res->p[i])
 672                         goto error;
 673                 res->n++;
 674         }
 675         for (; i < size; ++i) {
 676                 res->p[i] = isl_upoly_zero(up1->ctx);
 677                 if (!res->p[i])
 678                         goto error;
 679                 res->n++;
 680         }
 681         for (i = 0; i < rec1->n; ++i) {
 682                 for (j = 1; j < rec2->n; ++j) {
 683                         struct isl_upoly *up;
 684                         up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
 685                                             isl_upoly_copy(rec1->p[i]));
 686                         res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
 687                         if (!res->p[i + j])
 688                                 goto error;
 689                 }
 690         }
 691
 692         isl_upoly_free(up1);
 693         isl_upoly_free(up2);
 694
 695         return &res->up;
 696 error:
 697         isl_upoly_free(up1);
 698         isl_upoly_free(up2);
 699         isl_upoly_free(&res->up);
 700         return NULL;
 701 }
 702
 703 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
 704         __isl_take struct isl_upoly *up2)
 705 {
 706         if (!up1 || !up2)
 707                 goto error;
 708
 709         if (isl_upoly_is_nan(up1)) {
 710                 isl_upoly_free(up2);
 711                 return up1;
 712         }
 713
 714         if (isl_upoly_is_nan(up2)) {
 715                 isl_upoly_free(up1);
 716                 return up2;
 717         }
 718
 719         if (isl_upoly_is_zero(up1)) {
 720                 isl_upoly_free(up2);
 721                 return up1;
 722         }
 723
 724         if (isl_upoly_is_zero(up2)) {
 725                 isl_upoly_free(up1);
 726                 return up2;
 727         }
 728
 729         if (isl_upoly_is_one(up1)) {
 730                 isl_upoly_free(up1);
 731                 return up2;
 732         }
 733
 734         if (isl_upoly_is_one(up2)) {
 735                 isl_upoly_free(up2);
 736                 return up1;
 737         }
 738
 739         if (up1->var < up2->var)
 740                 return isl_upoly_mul(up2, up1);
 741
 742         if (up2->var < up1->var) {
 743                 int i;
 744                 struct isl_upoly_rec *rec;
 745                 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
 746                         isl_ctx *ctx = up1->ctx;
 747                         isl_upoly_free(up1);
 748                         isl_upoly_free(up2);
 749                         return isl_upoly_nan(ctx);
 750                 }
 751                 up1 = isl_upoly_cow(up1);
 752                 rec = isl_upoly_as_rec(up1);
 753                 if (!rec)
 754                         goto error;
 755
 756                 for (i = 0; i < rec->n; ++i) {
 757                         rec->p[i] = isl_upoly_mul(rec->p[i],
 758                                                     isl_upoly_copy(up2));
 759                         if (!rec->p[i])
 760                                 goto error;
 761                 }
 762                 isl_upoly_free(up2);
 763                 return up1;
 764         }
 765
 766         if (isl_upoly_is_cst(up1))
 767                 return isl_upoly_mul_cst(up1, up2);
 768
 769         return isl_upoly_mul_rec(up1, up2);
 770 error:
 771         isl_upoly_free(up1);
 772         isl_upoly_free(up2);
 773         return NULL;
 774 }
 775
 776 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
 777         unsigned n_div)
 778 {
 779         struct isl_qpolynomial *qp;
 780         unsigned total;
 781
 782         if (!dim)
 783                 return NULL;
 784
 785         total = isl_dim_total(dim);
 786
 787         qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
 788         if (!qp)
 789                 return NULL;
 790
 791         qp->ref = 1;
 792         qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
 793         if (!qp->div)
 794                 goto error;
 795
 796         qp->dim = dim;
 797
 798         return qp;
 799 error:
 800         isl_dim_free(dim);
 801         isl_qpolynomial_free(qp);
 802         return NULL;
 803 }
 804
 805 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
 806 {
 807         if (!qp)
 808                 return NULL;
 809
 810         qp->ref++;
 811         return qp;
 812 }
 813
 814 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
 815 {
 816         struct isl_qpolynomial *dup;
 817
 818         if (!qp)
 819                 return NULL;
 820
 821         dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row);
 822         if (!dup)
 823                 return NULL;
 824         isl_mat_free(dup->div);
 825         dup->div = isl_mat_copy(qp->div);
 826         if (!dup->div)
 827                 goto error;
 828         dup->upoly = isl_upoly_copy(qp->upoly);
 829         if (!dup->upoly)
 830                 goto error;
 831
 832         return dup;
 833 error:
 834         isl_qpolynomial_free(dup);
 835         return NULL;
 836 }
 837
 838 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
 839 {
 840         if (!qp)
 841                 return NULL;
 842
 843         if (qp->ref == 1)
 844                 return qp;
 845         qp->ref--;
 846         return isl_qpolynomial_dup(qp);
 847 }
 848
 849 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
 850 {
 851         if (!qp)
 852                 return;
 853
 854         if (--qp->ref > 0)
 855                 return;
 856
 857         isl_dim_free(qp->dim);
 858         isl_mat_free(qp->div);
 859         isl_upoly_free(qp->upoly);
 860
 861         free(qp);
 862 }
 863
 864 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
 865 {
 866         int n_row, n_col;
 867         int equal;
 868
 869         isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
 870                                 div1->n_col >= div2->n_col, return -1);
 871
 872         if (div1->n_row == div2->n_row)
 873                 return isl_mat_is_equal(div1, div2);
 874
 875         n_row = div1->n_row;
 876         n_col = div1->n_col;
 877         div1->n_row = div2->n_row;
 878         div1->n_col = div2->n_col;
 879
 880         equal = isl_mat_is_equal(div1, div2);
 881
 882         div1->n_row = n_row;
 883         div1->n_col = n_col;
 884
 885         return equal;
 886 }
 887
 888 static void expand_row(__isl_keep isl_mat *dst, int d,
 889         __isl_keep isl_mat *src, int s, int *exp)
 890 {
 891         int i;
 892         unsigned c = src->n_col - src->n_row;
 893
 894         isl_seq_cpy(dst->row[d], src->row[s], c);
 895         isl_seq_clr(dst->row[d] + c, dst->n_col - c);
 896
 897         for (i = 0; i < s; ++i)
 898                 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
 899 }
 900
 901 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
 902 {
 903         int li, lj;
 904
 905         li = isl_seq_last_non_zero(div->row[i], div->n_col);
 906         lj = isl_seq_last_non_zero(div->row[j], div->n_col);
 907
 908         if (li != lj)
 909                 return li - lj;
 910
 911         return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
 912 }
 913
 914 struct isl_div_sort_info {
 915         isl_mat *div;
 916         int      row;
 917 };
 918
 919 static int div_sort_cmp(const void *p1, const void *p2)
 920 {
 921         const struct isl_div_sort_info *i1, *i2;
 922         i1 = (const struct isl_div_sort_info *) p1;
 923         i2 = (const struct isl_div_sort_info *) p2;
 924
 925         return cmp_row(i1->div, i1->row, i2->row);
 926 }
 927
 928 static __isl_give isl_mat *sort_divs(__isl_take isl_mat *div)
 929 {
 930         int i;
 931         struct isl_div_sort_info *array = NULL;
 932         int *pos = NULL;
 933
 934         if (!div)
 935                 return NULL;
 936         if (div->n_row <= 1)
 937                 return div;
 938
 939         array = isl_alloc_array(div->ctx, struct isl_div_sort_info, div->n_row);
 940         pos = isl_alloc_array(div->ctx, int, div->n_row);
 941         if (!array || !pos)
 942                 goto error;
 943
 944         for (i = 0; i < div->n_row; ++i) {
 945                 array[i].div = div;
 946                 array[i].row = i;
 947                 pos[i] = i;
 948         }
 949
 950         qsort(array, div->n_row, sizeof(struct isl_div_sort_info),
 951                 div_sort_cmp);
 952
 953         for (i = 0; i < div->n_row; ++i) {
 954                 int t;
 955                 if (pos[array[i].row] == i)
 956                         continue;
 957                 div = isl_mat_cow(div);
 958                 div = isl_mat_swap_rows(div, i, pos[array[i].row]);
 959                 t = pos[array[i].row];
 960                 pos[array[i].row] = pos[i];
 961                 pos[i] = t;
 962         }
 963
 964         free(array);
 965
 966         return div;
 967 error:
 968         free(pos);
 969         free(array);
 970         isl_mat_free(div);
 971         return NULL;
 972 }
 973
 974 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
 975         __isl_keep isl_mat *div2, int *exp1, int *exp2)
 976 {
 977         int i, j, k;
 978         isl_mat *div = NULL;
 979         unsigned d = div1->n_col - div1->n_row;
 980
 981         div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
 982                                 d + div1->n_row + div2->n_row);
 983         if (!div)
 984                 return NULL;
 985
 986         for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
 987                 int cmp;
 988
 989                 expand_row(div, k, div1, i, exp1);
 990                 expand_row(div, k + 1, div2, j, exp2);
 991
 992                 cmp = cmp_row(div, k, k + 1);
 993                 if (cmp == 0) {
 994                         exp1[i++] = k;
 995                         exp2[j++] = k;
 996                 } else if (cmp < 0) {
 997                         exp1[i++] = k;
 998                 } else {
 999                         exp2[j++] = k;
1000                         isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
1001                 }
1002         }
1003         for (; i < div1->n_row; ++i, ++k) {
1004                 expand_row(div, k, div1, i, exp1);
1005                 exp1[i] = k;
1006         }
1007         for (; j < div2->n_row; ++j, ++k) {
1008                 expand_row(div, k, div2, j, exp2);
1009                 exp2[j] = k;
1010         }
1011
1012         div->n_row = k;
1013         div->n_col = d + k;
1014
1015         return div;
1016 }
1017
1018 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1019         int *exp, int first)
1020 {
1021         int i;
1022         struct isl_upoly_rec *rec;
1023
1024         if (isl_upoly_is_cst(up))
1025                 return up;
1026
1027         if (up->var < first)
1028                 return up;
1029
1030         if (exp[up->var - first] == up->var - first)
1031                 return up;
1032
1033         up = isl_upoly_cow(up);
1034         if (!up)
1035                 goto error;
1036
1037         up->var = exp[up->var - first] + first;
1038
1039         rec = isl_upoly_as_rec(up);
1040         if (!rec)
1041                 goto error;
1042
1043         for (i = 0; i < rec->n; ++i) {
1044                 rec->p[i] = expand(rec->p[i], exp, first);
1045                 if (!rec->p[i])
1046                         goto error;
1047         }
1048
1049         return up;
1050 error:
1051         isl_upoly_free(up);
1052         return NULL;
1053 }
1054
1055 static __isl_give isl_qpolynomial *with_merged_divs(
1056         __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1057                                           __isl_take isl_qpolynomial *qp2),
1058         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1059 {
1060         int *exp1 = NULL;
1061         int *exp2 = NULL;
1062         isl_mat *div = NULL;
1063
1064         qp1 = isl_qpolynomial_cow(qp1);
1065         qp2 = isl_qpolynomial_cow(qp2);
1066
1067         if (!qp1 || !qp2)
1068                 goto error;
1069
1070         isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1071                                 qp1->div->n_col >= qp2->div->n_col, goto error);
1072
1073         exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1074         exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1075         if (!exp1 || !exp2)
1076                 goto error;
1077
1078         div = merge_divs(qp1->div, qp2->div, exp1, exp2);
1079         if (!div)
1080                 goto error;
1081
1082         isl_mat_free(qp1->div);
1083         qp1->div = isl_mat_copy(div);
1084         isl_mat_free(qp2->div);
1085         qp2->div = isl_mat_copy(div);
1086
1087         qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1088         qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1089
1090         if (!qp1->upoly || !qp2->upoly)
1091                 goto error;
1092
1093         isl_mat_free(div);
1094         free(exp1);
1095         free(exp2);
1096
1097         return fn(qp1, qp2);
1098 error:
1099         isl_mat_free(div);
1100         free(exp1);
1101         free(exp2);
1102         isl_qpolynomial_free(qp1);
1103         isl_qpolynomial_free(qp2);
1104         return NULL;
1105 }
1106
1107 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1108         __isl_take isl_qpolynomial *qp2)
1109 {
1110         qp1 = isl_qpolynomial_cow(qp1);
1111
1112         if (!qp1 || !qp2)
1113                 goto error;
1114
1115         if (qp1->div->n_row < qp2->div->n_row)
1116                 return isl_qpolynomial_add(qp2, qp1);
1117
1118         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1119         if (!compatible_divs(qp1->div, qp2->div))
1120                 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1121
1122         qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1123         if (!qp1->upoly)
1124                 goto error;
1125
1126         isl_qpolynomial_free(qp2);
1127
1128         return qp1;
1129 error:
1130         isl_qpolynomial_free(qp1);
1131         isl_qpolynomial_free(qp2);
1132         return NULL;
1133 }
1134
1135 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1136         __isl_take isl_qpolynomial *qp2)
1137 {
1138         return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1139 }
1140
1141 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1142 {
1143         qp = isl_qpolynomial_cow(qp);
1144
1145         if (!qp)
1146                 return NULL;
1147
1148         qp->upoly = isl_upoly_neg(qp->upoly);
1149         if (!qp->upoly)
1150                 goto error;
1151
1152         return qp;
1153 error:
1154         isl_qpolynomial_free(qp);
1155         return NULL;
1156 }
1157
1158 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1159         __isl_take isl_qpolynomial *qp2)
1160 {
1161         qp1 = isl_qpolynomial_cow(qp1);
1162
1163         if (!qp1 || !qp2)
1164                 goto error;
1165
1166         if (qp1->div->n_row < qp2->div->n_row)
1167                 return isl_qpolynomial_mul(qp2, qp1);
1168
1169         isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1170         if (!compatible_divs(qp1->div, qp2->div))
1171                 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1172
1173         qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1174         if (!qp1->upoly)
1175                 goto error;
1176
1177         isl_qpolynomial_free(qp2);
1178
1179         return qp1;
1180 error:
1181         isl_qpolynomial_free(qp1);
1182         isl_qpolynomial_free(qp2);
1183         return NULL;
1184 }
1185
1186 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1187 {
1188         struct isl_qpolynomial *qp;
1189         struct isl_upoly_cst *cst;
1190
1191         qp = isl_qpolynomial_alloc(dim, 0);
1192         if (!qp)
1193                 return NULL;
1194
1195         qp->upoly = isl_upoly_zero(dim->ctx);
1196         if (!qp->upoly)
1197                 goto error;
1198
1199         return qp;
1200 error:
1201         isl_qpolynomial_free(qp);
1202         return NULL;
1203 }
1204
1205 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1206 {
1207         struct isl_qpolynomial *qp;
1208         struct isl_upoly_cst *cst;
1209
1210         qp = isl_qpolynomial_alloc(dim, 0);
1211         if (!qp)
1212                 return NULL;
1213
1214         qp->upoly = isl_upoly_infty(dim->ctx);
1215         if (!qp->upoly)
1216                 goto error;
1217
1218         return qp;
1219 error:
1220         isl_qpolynomial_free(qp);
1221         return NULL;
1222 }
1223
1224 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1225 {
1226         struct isl_qpolynomial *qp;
1227         struct isl_upoly_cst *cst;
1228
1229         qp = isl_qpolynomial_alloc(dim, 0);
1230         if (!qp)
1231                 return NULL;
1232
1233         qp->upoly = isl_upoly_nan(dim->ctx);
1234         if (!qp->upoly)
1235                 goto error;
1236
1237         return qp;
1238 error:
1239         isl_qpolynomial_free(qp);
1240         return NULL;
1241 }
1242
1243 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1244         isl_int v)
1245 {
1246         struct isl_qpolynomial *qp;
1247         struct isl_upoly_cst *cst;
1248
1249         qp = isl_qpolynomial_alloc(dim, 0);
1250         if (!qp)
1251                 return NULL;
1252
1253         qp->upoly = isl_upoly_zero(dim->ctx);
1254         if (!qp->upoly)
1255                 goto error;
1256         cst = isl_upoly_as_cst(qp->upoly);
1257         isl_int_set(cst->n, v);
1258
1259         return qp;
1260 error:
1261         isl_qpolynomial_free(qp);
1262         return NULL;
1263 }
1264
1265 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1266         isl_int *n, isl_int *d)
1267 {
1268         struct isl_upoly_cst *cst;
1269
1270         if (!qp)
1271                 return -1;
1272
1273         if (!isl_upoly_is_cst(qp->upoly))
1274                 return 0;
1275
1276         cst = isl_upoly_as_cst(qp->upoly);
1277         if (!cst)
1278                 return -1;
1279
1280         if (n)
1281                 isl_int_set(*n, cst->n);
1282         if (d)
1283                 isl_int_set(*d, cst->d);
1284
1285         return 1;
1286 }
1287
1288 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
1289         __isl_keep isl_qpolynomial *qp2)
1290 {
1291         if (!qp1 || !qp2)
1292                 return -1;
1293
1294         return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1295 }
1296
1297 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1298 {
1299         int i;
1300         struct isl_upoly_rec *rec;
1301
1302         if (isl_upoly_is_cst(up)) {
1303                 struct isl_upoly_cst *cst;
1304                 cst = isl_upoly_as_cst(up);
1305                 if (!cst)
1306                         return;
1307                 isl_int_lcm(*d, *d, cst->d);
1308                 return;
1309         }
1310
1311         rec = isl_upoly_as_rec(up);
1312         if (!rec)
1313                 return;
1314
1315         for (i = 0; i < rec->n; ++i)
1316                 upoly_update_den(rec->p[i], d);
1317 }
1318
1319 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1320 {
1321         isl_int_set_si(*d, 1);
1322         if (!qp)
1323                 return;
1324         upoly_update_den(qp->upoly, d);
1325 }
1326
1327 __isl_give struct isl_upoly *isl_upoly_pow(isl_ctx *ctx, int pos, int power)
1328 {
1329         int i;
1330         struct isl_upoly *up;
1331         struct isl_upoly_rec *rec;
1332         struct isl_upoly_cst *cst;
1333
1334         rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1335         if (!rec)
1336                 return NULL;
1337         for (i = 0; i < 1 + power; ++i) {
1338                 rec->p[i] = isl_upoly_zero(ctx);
1339                 if (!rec->p[i])
1340                         goto error;
1341                 rec->n++;
1342         }
1343         cst = isl_upoly_as_cst(rec->p[power]);
1344         isl_int_set_si(cst->n, 1);
1345
1346         return &rec->up;
1347 error:
1348         isl_upoly_free(&rec->up);
1349         return NULL;
1350 }
1351
1352 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1353         int pos, int power)
1354 {
1355         struct isl_qpolynomial *qp;
1356         struct isl_ctx *ctx;
1357
1358         if (!dim)
1359                 return NULL;
1360
1361         ctx = dim->ctx;
1362
1363         qp = isl_qpolynomial_alloc(dim, 0);
1364         if (!qp)
1365                 return NULL;
1366
1367         qp->upoly = isl_upoly_pow(ctx, pos, power);
1368         if (!qp->upoly)
1369                 goto error;
1370
1371         return qp;
1372 error:
1373         isl_qpolynomial_free(qp);
1374         return NULL;
1375 }
1376
1377 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1378         enum isl_dim_type type, unsigned pos)
1379 {
1380         if (!dim)
1381                 return NULL;
1382
1383         isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1384         isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1385
1386         if (type == isl_dim_set)
1387                 pos += isl_dim_size(dim, isl_dim_param);
1388
1389         return isl_qpolynomial_pow(dim, pos, 1);
1390 error:
1391         isl_dim_free(dim);
1392         return NULL;
1393 }
1394
1395 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1396         int power)
1397 {
1398         struct isl_qpolynomial *qp = NULL;
1399         struct isl_upoly_rec *rec;
1400         struct isl_upoly_cst *cst;
1401         int i;
1402         int pos;
1403
1404         if (!div)
1405                 return NULL;
1406         isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1407
1408         qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1);
1409         if (!qp)
1410                 goto error;
1411
1412         isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1413         isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1414
1415         pos = isl_dim_total(qp->dim);
1416         qp->upoly = &isl_upoly_alloc_rec(div->ctx, pos, 1 + power)->up;
1417         if (!qp->upoly)
1418                 goto error;
1419         rec = isl_upoly_as_rec(qp->upoly);
1420         for (i = 0; i < 1 + power; ++i) {
1421                 rec->p[i] = isl_upoly_zero(div->ctx);
1422                 if (!rec->p[i])
1423                         goto error;
1424                 rec->n++;
1425         }
1426         cst = isl_upoly_as_cst(rec->p[power]);
1427         isl_int_set_si(cst->n, 1);
1428
1429         isl_div_free(div);
1430
1431         return qp;
1432 error:
1433         isl_qpolynomial_free(qp);
1434         isl_div_free(div);
1435         return NULL;
1436 }
1437
1438 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1439 {
1440         return isl_qpolynomial_div_pow(div, 1);
1441 }
1442
1443 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1444         const isl_int n, const isl_int d)
1445 {
1446         struct isl_qpolynomial *qp;
1447         struct isl_upoly_cst *cst;
1448
1449         qp = isl_qpolynomial_alloc(dim, 0);
1450         if (!qp)
1451                 return NULL;
1452
1453         qp->upoly = isl_upoly_zero(dim->ctx);
1454         if (!qp->upoly)
1455                 goto error;
1456         cst = isl_upoly_as_cst(qp->upoly);
1457         isl_int_set(cst->n, n);
1458         isl_int_set(cst->d, d);
1459
1460         return qp;
1461 error:
1462         isl_qpolynomial_free(qp);
1463         return NULL;
1464 }
1465
1466 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
1467 {
1468         struct isl_upoly_rec *rec;
1469         int i;
1470
1471         if (!up)
1472                 return -1;
1473
1474         if (isl_upoly_is_cst(up))
1475                 return 0;
1476
1477         if (up->var < d)
1478                 active[up->var] = 1;
1479
1480         rec = isl_upoly_as_rec(up);
1481         for (i = 0; i < rec->n; ++i)
1482                 if (up_set_active(rec->p[i], active, d) < 0)
1483                         return -1;
1484
1485         return 0;
1486 }
1487
1488 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
1489 {
1490         int i, j;
1491         int d = isl_dim_total(qp->dim);
1492
1493         if (!qp || !active)
1494                 return -1;
1495
1496         for (i = 0; i < d; ++i)
1497                 for (j = 0; j < qp->div->n_row; ++j) {
1498                         if (isl_int_is_zero(qp->div->row[j][2 + i]))
1499                                 continue;
1500                         active[i] = 1;
1501                         break;
1502                 }
1503
1504         return up_set_active(qp->upoly, active, d);
1505 }
1506
1507 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
1508         enum isl_dim_type type, unsigned first, unsigned n)
1509 {
1510         int i;
1511         int *active = NULL;
1512         int involves = 0;
1513
1514         if (!qp)
1515                 return -1;
1516         if (n == 0)
1517                 return 0;
1518
1519         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1520                         return -1);
1521         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1522                                  type == isl_dim_set, return -1);
1523
1524         active = isl_calloc_array(set->ctx, int, isl_dim_total(qp->dim));
1525         if (set_active(qp, active) < 0)
1526                 goto error;
1527
1528         if (type == isl_dim_set)
1529                 first += isl_dim_size(qp->dim, isl_dim_param);
1530         for (i = 0; i < n; ++i)
1531                 if (active[first + i]) {
1532                         involves = 1;
1533                         break;
1534                 }
1535
1536         free(active);
1537
1538         return involves;
1539 error:
1540         free(active);
1541         return -1;
1542 }
1543
1544 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
1545         enum isl_dim_type type, unsigned first, unsigned n)
1546 {
1547         int i;
1548
1549         if (!fold)
1550                 return -1;
1551         if (fold->n == 0 || n == 0)
1552                 return 0;
1553
1554         for (i = 0; i < fold->n; ++i) {
1555                 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
1556                                                             type, first, n);
1557                 if (involves < 0 || involves)
1558                         return involves;
1559         }
1560         return 0;
1561 }
1562
1563 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
1564         unsigned first, unsigned n)
1565 {
1566         int i;
1567         struct isl_upoly_rec *rec;
1568
1569         if (!up)
1570                 return NULL;
1571         if (n == 0 || up->var < 0 || up->var < first)
1572                 return up;
1573         if (up->var < first + n) {
1574                 up = replace_by_constant_term(up);
1575                 return isl_upoly_drop(up, first, n);
1576         }
1577         up = isl_upoly_cow(up);
1578         if (!up)
1579                 return NULL;
1580         up->var -= n;
1581         rec = isl_upoly_as_rec(up);
1582         if (!rec)
1583                 goto error;
1584
1585         for (i = 0; i < rec->n; ++i) {
1586                 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
1587                 if (!rec->p[i])
1588                         goto error;
1589         }
1590
1591         return up;
1592 error:
1593         isl_upoly_free(up);
1594         return NULL;
1595 }
1596
1597 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
1598         __isl_take isl_qpolynomial *qp,
1599         enum isl_dim_type type, unsigned first, unsigned n)
1600 {
1601         if (!qp)
1602                 return NULL;
1603         if (n == 0)
1604                 return qp;
1605
1606         qp = isl_qpolynomial_cow(qp);
1607         if (!qp)
1608                 return NULL;
1609
1610         isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1611                         goto error);
1612         isl_assert(qp->dim->ctx, type == isl_dim_param ||
1613                                  type == isl_dim_set, goto error);
1614
1615         qp->dim = isl_dim_drop(qp->dim, type, first, n);
1616         if (!qp->dim)
1617                 goto error;
1618
1619         if (type == isl_dim_set)
1620                 first += isl_dim_size(qp->dim, isl_dim_param);
1621
1622         qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
1623         if (!qp->div)
1624                 goto error;
1625
1626         qp->upoly = isl_upoly_drop(qp->upoly, first, n);
1627         if (!qp->upoly)
1628                 goto error;
1629
1630         return qp;
1631 error:
1632         isl_qpolynomial_free(qp);
1633         return NULL;
1634 }
1635
1636 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
1637         __isl_take isl_qpolynomial_fold *fold,
1638         enum isl_dim_type type, unsigned first, unsigned n)
1639 {
1640         int i;
1641
1642         if (!fold)
1643                 return NULL;
1644         if (n == 0)
1645                 return fold;
1646
1647         fold = isl_qpolynomial_fold_cow(fold);
1648         if (!fold)
1649                 return NULL;
1650         fold->dim = isl_dim_drop(fold->dim, type, first, n);
1651         if (!fold->dim)
1652                 goto error;
1653
1654         for (i = 0; i < fold->n; ++i) {
1655                 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
1656                                                             type, first, n);
1657                 if (!fold->qp[i])
1658                         goto error;
1659         }
1660
1661         return fold;
1662 error:
1663         isl_qpolynomial_fold_free(fold);
1664         return NULL;
1665 }
1666
1667 #undef PW
1668 #define PW isl_pw_qpolynomial
1669 #undef EL
1670 #define EL isl_qpolynomial
1671 #undef IS_ZERO
1672 #define IS_ZERO is_zero
1673 #undef FIELD
1674 #define FIELD qp
1675 #undef ADD
1676 #define ADD add
1677
1678 #include <isl_pw_templ.c>
1679
1680 #undef PW
1681 #define PW isl_pw_qpolynomial_fold
1682 #undef EL
1683 #define EL isl_qpolynomial_fold
1684 #undef IS_ZERO
1685 #define IS_ZERO is_empty
1686 #undef FIELD
1687 #define FIELD fold
1688 #undef ADD
1689 #define ADD fold
1690
1691 #include <isl_pw_templ.c>
1692
1693 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1694 {
1695         if (!pwqp)
1696                 return -1;
1697
1698         if (pwqp->n != -1)
1699                 return 0;
1700
1701         if (!isl_set_fast_is_universe(pwqp->p[0].set))
1702                 return 0;
1703
1704         return isl_qpolynomial_is_one(pwqp->p[0].qp);
1705 }
1706
1707 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1708         __isl_take isl_pw_qpolynomial *pwqp1,
1709         __isl_take isl_pw_qpolynomial *pwqp2)
1710 {
1711         int i, j, n;
1712         struct isl_pw_qpolynomial *res;
1713         isl_set *set;
1714
1715         if (!pwqp1 || !pwqp2)
1716                 goto error;
1717
1718         isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1719                         goto error);
1720
1721         if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1722                 isl_pw_qpolynomial_free(pwqp2);
1723                 return pwqp1;
1724         }
1725
1726         if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1727                 isl_pw_qpolynomial_free(pwqp1);
1728                 return pwqp2;
1729         }
1730
1731         if (isl_pw_qpolynomial_is_one(pwqp1)) {
1732                 isl_pw_qpolynomial_free(pwqp1);
1733                 return pwqp2;
1734         }
1735
1736         if (isl_pw_qpolynomial_is_one(pwqp2)) {
1737                 isl_pw_qpolynomial_free(pwqp2);
1738                 return pwqp1;
1739         }
1740
1741         n = pwqp1->n * pwqp2->n;
1742         res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
1743
1744         for (i = 0; i < pwqp1->n; ++i) {
1745                 for (j = 0; j < pwqp2->n; ++j) {
1746                         struct isl_set *common;
1747                         struct isl_qpolynomial *prod;
1748                         common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1749                                                 isl_set_copy(pwqp2->p[j].set));
1750                         if (isl_set_fast_is_empty(common)) {
1751                                 isl_set_free(common);
1752                                 continue;
1753                         }
1754
1755                         prod = isl_qpolynomial_mul(
1756                                 isl_qpolynomial_copy(pwqp1->p[i].qp),
1757                                 isl_qpolynomial_copy(pwqp2->p[j].qp));
1758
1759                         res = isl_pw_qpolynomial_add_piece(res, common, prod);
1760                 }
1761         }
1762
1763         isl_pw_qpolynomial_free(pwqp1);
1764         isl_pw_qpolynomial_free(pwqp2);
1765
1766         return res;
1767 error:
1768         isl_pw_qpolynomial_free(pwqp1);
1769         isl_pw_qpolynomial_free(pwqp2);
1770         return NULL;
1771 }
1772
1773 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1774         __isl_take isl_pw_qpolynomial *pwqp)
1775 {
1776         int i, j, n;
1777         struct isl_pw_qpolynomial *res;
1778         isl_set *set;
1779
1780         if (!pwqp)
1781                 return NULL;
1782
1783         if (isl_pw_qpolynomial_is_zero(pwqp))
1784                 return pwqp;
1785
1786         pwqp = isl_pw_qpolynomial_cow(pwqp);
1787
1788         for (i = 0; i < pwqp->n; ++i) {
1789                 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1790                 if (!pwqp->p[i].qp)
1791                         goto error;
1792         }
1793
1794         return pwqp;
1795 error:
1796         isl_pw_qpolynomial_free(pwqp);
1797         return NULL;
1798 }
1799
1800 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1801         __isl_take isl_pw_qpolynomial *pwqp1,
1802         __isl_take isl_pw_qpolynomial *pwqp2)
1803 {
1804         return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1805 }
1806
1807 __isl_give struct isl_upoly *isl_upoly_eval(
1808         __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1809 {
1810         int i;
1811         struct isl_upoly_rec *rec;
1812         struct isl_upoly *res;
1813         struct isl_upoly *base;
1814
1815         if (isl_upoly_is_cst(up)) {
1816                 isl_vec_free(vec);
1817                 return up;
1818         }
1819
1820         rec = isl_upoly_as_rec(up);
1821         if (!rec)
1822                 goto error;
1823
1824         isl_assert(up->ctx, rec->n >= 1, goto error);
1825
1826         base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1827
1828         res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1829                                 isl_vec_copy(vec));
1830
1831         for (i = rec->n - 2; i >= 0; --i) {
1832                 res = isl_upoly_mul(res, isl_upoly_copy(base));
1833                 res = isl_upoly_sum(res,
1834                             isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1835                                                             isl_vec_copy(vec)));
1836         }
1837
1838         isl_upoly_free(base);
1839         isl_upoly_free(up);
1840         isl_vec_free(vec);
1841         return res;
1842 error:
1843         isl_upoly_free(up);
1844         isl_vec_free(vec);
1845         return NULL;
1846 }
1847
1848 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1849         __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1850 {
1851         isl_vec *ext;
1852         struct isl_upoly *up;
1853         isl_dim *dim;
1854
1855         if (!qp || !pnt)
1856                 goto error;
1857         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1858
1859         if (qp->div->n_row == 0)
1860                 ext = isl_vec_copy(pnt->vec);
1861         else {
1862                 int i;
1863                 unsigned dim = isl_dim_total(qp->dim);
1864                 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1865                 if (!ext)
1866                         goto error;
1867
1868                 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1869                 for (i = 0; i < qp->div->n_row; ++i) {
1870                         isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1871                                                 1 + dim + i, &ext->el[1+dim+i]);
1872                         isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1873                                         qp->div->row[i][0]);
1874                 }
1875         }
1876
1877         up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1878         if (!up)
1879                 goto error;
1880
1881         dim = isl_dim_copy(qp->dim);
1882         isl_qpolynomial_free(qp);
1883         isl_point_free(pnt);
1884
1885         qp = isl_qpolynomial_alloc(dim, 0);
1886         if (!qp)
1887                 isl_upoly_free(up);
1888         else
1889                 qp->upoly = up;
1890
1891         return qp;
1892 error:
1893         isl_qpolynomial_free(qp);
1894         isl_point_free(pnt);
1895         return NULL;
1896 }
1897
1898 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
1899         __isl_keep struct isl_upoly_cst *cst2)
1900 {
1901         int cmp;
1902         isl_int t;
1903         isl_int_init(t);
1904         isl_int_mul(t, cst1->n, cst2->d);
1905         isl_int_submul(t, cst2->n, cst1->d);
1906         cmp = isl_int_sgn(t);
1907         isl_int_clear(t);
1908         return cmp;
1909 }
1910
1911 static __isl_give isl_qpolynomial *qpolynomial_min(
1912         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1913 {
1914         struct isl_upoly_cst *cst1, *cst2;
1915         int cmp;
1916
1917         if (!qp1 || !qp2)
1918                 goto error;
1919         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1920         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1921         cst1 = isl_upoly_as_cst(qp1->upoly);
1922         cst2 = isl_upoly_as_cst(qp2->upoly);
1923         cmp = isl_upoly_cmp(cst1, cst2);
1924
1925         if (cmp <= 0) {
1926                 isl_qpolynomial_free(qp2);
1927         } else {
1928                 isl_qpolynomial_free(qp1);
1929                 qp1 = qp2;
1930         }
1931         return qp1;
1932 error:
1933         isl_qpolynomial_free(qp1);
1934         isl_qpolynomial_free(qp2);
1935         return NULL;
1936 }
1937
1938 static __isl_give isl_qpolynomial *qpolynomial_max(
1939         __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1940 {
1941         struct isl_upoly_cst *cst1, *cst2;
1942         int cmp;
1943
1944         if (!qp1 || !qp2)
1945                 goto error;
1946         isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1947         isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1948         cst1 = isl_upoly_as_cst(qp1->upoly);
1949         cst2 = isl_upoly_as_cst(qp2->upoly);
1950         cmp = isl_upoly_cmp(cst1, cst2);
1951
1952         if (cmp >= 0) {
1953                 isl_qpolynomial_free(qp2);
1954         } else {
1955                 isl_qpolynomial_free(qp1);
1956                 qp1 = qp2;
1957         }
1958         return qp1;
1959 error:
1960         isl_qpolynomial_free(qp1);
1961         isl_qpolynomial_free(qp2);
1962         return NULL;
1963 }
1964
1965 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
1966         __isl_keep isl_qpolynomial_fold *fold2)
1967 {
1968         int i;
1969
1970         if (!fold1 || !fold2)
1971                 return -1;
1972
1973         if (fold1->n != fold2->n)
1974                 return 0;
1975
1976         /* We probably want to sort the qps first... */
1977         for (i = 0; i < fold1->n; ++i) {
1978                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
1979                 if (eq < 0 || !eq)
1980                         return eq;
1981         }
1982
1983         return 1;
1984 }
1985
1986 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
1987         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
1988 {
1989         isl_qpolynomial *qp;
1990
1991         if (!fold || !pnt)
1992                 goto error;
1993         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
1994         isl_assert(pnt->dim->ctx,
1995                 fold->type == isl_fold_max || fold->type == isl_fold_min,
1996                 goto error);
1997
1998         if (fold->n == 0)
1999                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
2000         else {
2001                 int i;
2002                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
2003                                                 isl_point_copy(pnt));
2004                 for (i = 1; i < fold->n; ++i) {
2005                         isl_qpolynomial *qp_i;
2006                         qp_i = isl_qpolynomial_eval(
2007                                             isl_qpolynomial_copy(fold->qp[i]),
2008                                             isl_point_copy(pnt));
2009                         if (fold->type == isl_fold_max)
2010                                 qp = qpolynomial_max(qp, qp_i);
2011                         else
2012                                 qp = qpolynomial_min(qp, qp_i);
2013                 }
2014         }
2015         isl_qpolynomial_fold_free(fold);
2016         isl_point_free(pnt);
2017
2018         return qp;
2019 error:
2020         isl_qpolynomial_fold_free(fold);
2021         isl_point_free(pnt);
2022         return NULL;
2023 }
2024
2025 static int *reordering_move(isl_ctx *ctx,
2026         unsigned len, unsigned dst, unsigned src, unsigned n)
2027 {
2028         int i;
2029         int *reordering;
2030
2031         reordering = isl_alloc_array(ctx, int, len);
2032         if (!reordering)
2033                 return NULL;
2034
2035         if (dst <= src) {
2036                 for (i = 0; i < dst; ++i)
2037                         reordering[i] = i;
2038                 for (i = 0; i < n; ++i)
2039                         reordering[src + i] = dst + i;
2040                 for (i = 0; i < src - dst; ++i)
2041                         reordering[dst + i] = dst + n + i;
2042                 for (i = 0; i < len - src - n; ++i)
2043                         reordering[src + n + i] = src + n + i;
2044         } else {
2045                 for (i = 0; i < src; ++i)
2046                         reordering[i] = i;
2047                 for (i = 0; i < n; ++i)
2048                         reordering[src + i] = dst + i;
2049                 for (i = 0; i < dst - src; ++i)
2050                         reordering[src + n + i] = src + i;
2051                 for (i = 0; i < len - dst - n; ++i)
2052                         reordering[dst + n + i] = dst + n + i;
2053         }
2054
2055         return reordering;
2056 }
2057
2058 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
2059         int *r)
2060 {
2061         int i;
2062         struct isl_upoly_rec *rec;
2063         struct isl_upoly *base;
2064         struct isl_upoly *res;
2065
2066         if (isl_upoly_is_cst(up))
2067                 return up;
2068
2069         rec = isl_upoly_as_rec(up);
2070         if (!rec)
2071                 goto error;
2072
2073         isl_assert(up->ctx, rec->n >= 1, goto error);
2074
2075         base = isl_upoly_pow(up->ctx, r[up->var], 1);
2076         res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
2077
2078         for (i = rec->n - 2; i >= 0; --i) {
2079                 res = isl_upoly_mul(res, isl_upoly_copy(base));
2080                 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
2081         }
2082
2083         isl_upoly_free(base);
2084         isl_upoly_free(up);
2085
2086         return res;
2087 error:
2088         isl_upoly_free(up);
2089         return NULL;
2090 }
2091
2092 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
2093         enum isl_dim_type dst_type, unsigned dst_pos,
2094         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
2095 {
2096         unsigned g_dst_pos;
2097         unsigned g_src_pos;
2098         int *reordering;
2099
2100         qp = isl_qpolynomial_cow(qp);
2101         if (!qp)
2102                 return NULL;
2103
2104         isl_assert(qp->dim->ctx, src_pos + n <= isl_dim_size(qp->dim, src_type),
2105                 goto error);
2106
2107         g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
2108         g_src_pos = pos(qp->dim, src_type) + src_pos;
2109         if (dst_type > src_type)
2110                 g_dst_pos -= n;
2111
2112         qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
2113         qp->div = sort_divs(qp->div);
2114         if (!qp->div)
2115                 goto error;
2116
2117         reordering = reordering_move(qp->dim->ctx,
2118                                 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
2119         if (!reordering)
2120                 goto error;
2121
2122         qp->upoly = reorder(qp->upoly, reordering);
2123         free(reordering);
2124         if (!qp->upoly)
2125                 goto error;
2126
2127         qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
2128         if (!qp->dim)
2129                 goto error;
2130
2131         return qp;
2132 error:
2133         isl_qpolynomial_free(qp);
2134         return NULL;
2135 }
2136
2137 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
2138         __isl_take isl_pw_qpolynomial *pwqp,
2139         enum isl_dim_type dst_type, unsigned dst_pos,
2140         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
2141 {
2142         int i;
2143
2144         pwqp = isl_pw_qpolynomial_cow(pwqp);
2145         if (!pwqp)
2146                 return NULL;
2147
2148         pwqp->dim = isl_dim_move(pwqp->dim,
2149                                     dst_type, dst_pos, src_type, src_pos, n);
2150         if (!pwqp->dim)
2151                 goto error;
2152
2153         for (i = 0; i < pwqp->n; ++i) {
2154                 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
2155                                                 src_type, src_pos, n);
2156                 if (!pwqp->p[i].set)
2157                         goto error;
2158                 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
2159                                         dst_type, dst_pos, src_type, src_pos, n);
2160                 if (!pwqp->p[i].qp)
2161                         goto error;
2162         }
2163
2164         return pwqp;
2165 error:
2166         isl_pw_qpolynomial_free(pwqp);
2167         return NULL;
2168 }
2169
2170 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
2171         __isl_keep isl_pw_qpolynomial *pwqp)
2172 {
2173         if (!pwqp)
2174                 return NULL;
2175
2176         return isl_dim_copy(pwqp->dim);
2177 }
2178
2179 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
2180         enum isl_dim_type type)
2181 {
2182         return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
2183 }
2184
2185 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
2186         __isl_take isl_mat *div)
2187 {
2188         isl_term *term;
2189         int n;
2190
2191         if (!dim || !div)
2192                 goto error;
2193
2194         n = isl_dim_total(dim) + div->n_row;
2195
2196         term = isl_calloc(dim->ctx, struct isl_term,
2197                         sizeof(struct isl_term) + (n - 1) * sizeof(int));
2198         if (!term)
2199                 goto error;
2200
2201         term->ref = 1;
2202         term->dim = dim;
2203         term->div = div;
2204         isl_int_init(term->n);
2205         isl_int_init(term->d);
2206
2207         return term;
2208 error:
2209         isl_dim_free(dim);
2210         isl_mat_free(div);
2211         return NULL;
2212 }
2213
2214 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
2215 {
2216         if (!term)
2217                 return NULL;
2218
2219         term->ref++;
2220         return term;
2221 }
2222
2223 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
2224 {
2225         int i;
2226         isl_term *dup;
2227         unsigned total;
2228
2229         if (term)
2230                 return NULL;
2231
2232         total = isl_dim_total(term->dim) + term->div->n_row;
2233
2234         dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
2235         if (!dup)
2236                 return NULL;
2237
2238         isl_int_set(dup->n, term->n);
2239         isl_int_set(dup->d, term->d);
2240
2241         for (i = 0; i < total; ++i)
2242                 dup->pow[i] = term->pow[i];
2243
2244         return dup;
2245 }
2246
2247 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
2248 {
2249         if (!term)
2250                 return NULL;
2251
2252         if (term->ref == 1)
2253                 return term;
2254         term->ref--;
2255         return isl_term_dup(term);
2256 }
2257
2258 void isl_term_free(__isl_take isl_term *term)
2259 {
2260         if (!term)
2261                 return;
2262
2263         if (--term->ref > 0)
2264                 return;
2265
2266         isl_dim_free(term->dim);
2267         isl_mat_free(term->div);
2268         isl_int_clear(term->n);
2269         isl_int_clear(term->d);
2270         free(term);
2271 }
2272
2273 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
2274 {
2275         if (!term)
2276                 return 0;
2277
2278         switch (type) {
2279         case isl_dim_param:
2280         case isl_dim_in:
2281         case isl_dim_out:       return isl_dim_size(term->dim, type);
2282         case isl_dim_div:       return term->div->n_row;
2283         case isl_dim_all:       return isl_dim_total(term->dim) + term->div->n_row;
2284         }
2285 }
2286
2287 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
2288 {
2289         return term ? term->dim->ctx : NULL;
2290 }
2291
2292 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
2293 {
2294         if (!term)
2295                 return;
2296         isl_int_set(*n, term->n);
2297 }
2298
2299 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
2300 {
2301         if (!term)
2302                 return;
2303         isl_int_set(*d, term->d);
2304 }
2305
2306 int isl_term_get_exp(__isl_keep isl_term *term,
2307         enum isl_dim_type type, unsigned pos)
2308 {
2309         if (!term)
2310                 return -1;
2311
2312         isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
2313
2314         if (type >= isl_dim_set)
2315                 pos += isl_dim_size(term->dim, isl_dim_param);
2316         if (type >= isl_dim_div)
2317                 pos += isl_dim_size(term->dim, isl_dim_set);
2318
2319         return term->pow[pos];
2320 }
2321
2322 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
2323 {
2324         isl_basic_map *bmap;
2325         unsigned total;
2326         int k;
2327
2328         if (!term)
2329                 return NULL;
2330
2331         isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
2332                         return NULL);
2333
2334         total = term->div->n_col - term->div->n_row - 2;
2335         /* No nested divs for now */
2336         isl_assert(term->dim->ctx,
2337                 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
2338                                         term->div->n_row) == -1,
2339                 return NULL);
2340
2341         bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
2342         if ((k = isl_basic_map_alloc_div(bmap)) < 0)
2343                 goto error;
2344
2345         isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
2346
2347         return isl_basic_map_div(bmap, k);
2348 error:
2349         isl_basic_map_free(bmap);
2350         return NULL;
2351 }
2352
2353 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
2354         int (*fn)(__isl_take isl_term *term, void *user),
2355         __isl_take isl_term *term, void *user)
2356 {
2357         int i;
2358         struct isl_upoly_rec *rec;
2359
2360         if (!up || !term)
2361                 goto error;
2362
2363         if (isl_upoly_is_zero(up))
2364                 return term;
2365
2366         isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
2367         isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
2368         isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
2369
2370         if (isl_upoly_is_cst(up)) {
2371                 struct isl_upoly_cst *cst;
2372                 cst = isl_upoly_as_cst(up);
2373                 if (!cst)
2374                         goto error;
2375                 term = isl_term_cow(term);
2376                 if (!term)
2377                         goto error;
2378                 isl_int_set(term->n, cst->n);
2379                 isl_int_set(term->d, cst->d);
2380                 if (fn(isl_term_copy(term), user) < 0)
2381                         goto error;
2382                 return term;
2383         }
2384
2385         rec = isl_upoly_as_rec(up);
2386         if (!rec)
2387                 goto error;
2388
2389         for (i = 0; i < rec->n; ++i) {
2390                 term = isl_term_cow(term);
2391                 if (!term)
2392                         goto error;
2393                 term->pow[up->var] = i;
2394                 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2395                 if (!term)
2396                         goto error;
2397         }
2398         term->pow[up->var] = 0;
2399
2400         return term;
2401 error:
2402         isl_term_free(term);
2403         return NULL;
2404 }
2405
2406 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
2407         int (*fn)(__isl_take isl_term *term, void *user), void *user)
2408 {
2409         isl_term *term;
2410
2411         if (!qp)
2412                 return -1;
2413
2414         term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
2415         if (!term)
2416                 return -1;
2417
2418         term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
2419
2420         isl_term_free(term);
2421
2422         return term ? 0 : -1;
2423 }
2424
2425 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2426         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2427                     void *user), void *user)
2428 {
2429         int i;
2430
2431         if (!pwqp)
2432                 return -1;
2433
2434         for (i = 0; i < pwqp->n; ++i)
2435                 if (fn(isl_set_copy(pwqp->p[i].set),
2436                                 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
2437                         return -1;
2438
2439         return 0;
2440 }
2441
2442 static int any_divs(__isl_keep isl_set *set)
2443 {
2444         int i;
2445
2446         if (!set)
2447                 return -1;
2448
2449         for (i = 0; i < set->n; ++i)
2450                 if (set->p[i]->n_div > 0)
2451                         return 1;
2452
2453         return 0;
2454 }
2455
2456 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
2457         __isl_take isl_dim *dim)
2458 {
2459         if (!qp || !dim)
2460                 goto error;
2461
2462         if (isl_dim_equal(qp->dim, dim)) {
2463                 isl_dim_free(dim);
2464                 return qp;
2465         }
2466
2467         qp = isl_qpolynomial_cow(qp);
2468         if (!qp)
2469                 goto error;
2470
2471         if (qp->div->n_row) {
2472                 int i;
2473                 int extra;
2474                 unsigned total;
2475                 int *exp;
2476
2477                 extra = isl_dim_size(dim, isl_dim_set) -
2478                                 isl_dim_size(qp->dim, isl_dim_set);
2479                 total = isl_dim_total(qp->dim);
2480                 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
2481                 if (!exp)
2482                         goto error;
2483                 for (i = 0; i < qp->div->n_row; ++i)
2484                         exp[i] = extra + i;
2485                 qp->upoly = expand(qp->upoly, exp, total);
2486                 free(exp);
2487                 if (!qp->upoly)
2488                         goto error;
2489                 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
2490                 if (!qp->div)
2491                         goto error;
2492                 for (i = 0; i < qp->div->n_row; ++i)
2493                         isl_seq_clr(qp->div->row[i] + 2 + total, extra);
2494         }
2495
2496         isl_dim_free(qp->dim);
2497         qp->dim = dim;
2498
2499         return qp;
2500 error:
2501         isl_dim_free(dim);
2502         isl_qpolynomial_free(qp);
2503         return NULL;
2504 }
2505
2506 static int foreach_lifted_subset(__isl_take isl_set *set,
2507         __isl_take isl_qpolynomial *qp,
2508         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2509                     void *user), void *user)
2510 {
2511         int i;
2512
2513         if (!set || !qp)
2514                 goto error;
2515
2516         for (i = 0; i < set->n; ++i) {
2517                 isl_set *lift;
2518                 isl_qpolynomial *copy;
2519
2520                 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2521                 lift = isl_set_lift(lift);
2522
2523                 copy = isl_qpolynomial_copy(qp);
2524                 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2525
2526                 if (fn(lift, copy, user) < 0)
2527                         goto error;
2528         }
2529
2530         isl_set_free(set);
2531         isl_qpolynomial_free(qp);
2532
2533         return 0;
2534 error:
2535         isl_set_free(set);
2536         isl_qpolynomial_free(qp);
2537         return -1;
2538 }
2539
2540 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2541         int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2542                     void *user), void *user)
2543 {
2544         int i;
2545
2546         if (!pwqp)
2547                 return -1;
2548
2549         for (i = 0; i < pwqp->n; ++i) {
2550                 isl_set *set;
2551                 isl_qpolynomial *qp;
2552
2553                 set = isl_set_copy(pwqp->p[i].set);
2554                 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2555                 if (!any_divs(set)) {
2556                         if (fn(set, qp, user) < 0)
2557                                 return -1;
2558                         continue;
2559                 }
2560                 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2561                         return -1;
2562         }
2563
2564         return 0;
2565 }
2566
2567 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
2568         enum isl_fold type, __isl_take isl_dim *dim, int n)
2569 {
2570         isl_qpolynomial_fold *fold;
2571
2572         if (!dim)
2573                 goto error;
2574
2575         isl_assert(dim->ctx, n >= 0, goto error);
2576         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
2577                         sizeof(struct isl_qpolynomial_fold) +
2578                         (n - 1) * sizeof(struct isl_qpolynomial *));
2579         if (!fold)
2580                 goto error;
2581
2582         fold->ref = 1;
2583         fold->size = n;
2584         fold->n = 0;
2585         fold->type = type;
2586         fold->dim = dim;
2587
2588         return fold;
2589 error:
2590         isl_dim_free(dim);
2591         return NULL;
2592 }
2593
2594 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
2595         __isl_take isl_dim *dim)
2596 {
2597         return qpolynomial_fold_alloc(type, dim, 0);
2598 }
2599
2600 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
2601         enum isl_fold type, __isl_take isl_qpolynomial *qp)
2602 {
2603         isl_qpolynomial_fold *fold;
2604
2605         if (!qp)
2606                 return NULL;
2607
2608         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
2609         if (!fold)
2610                 goto error;
2611
2612         fold->qp[0] = qp;
2613         fold->n++;
2614
2615         return fold;
2616 error:
2617         isl_qpolynomial_fold_free(fold);
2618         isl_qpolynomial_free(qp);
2619         return NULL;
2620 }
2621
2622 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
2623         __isl_keep isl_qpolynomial_fold *fold)
2624 {
2625         if (!fold)
2626                 return NULL;
2627
2628         fold->ref++;
2629         return fold;
2630 }
2631
2632 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
2633         __isl_keep isl_qpolynomial_fold *fold)
2634 {
2635         int i;
2636         isl_qpolynomial_fold *dup;
2637
2638         if (!fold)
2639                 return NULL;
2640         dup = qpolynomial_fold_alloc(fold->type,
2641                                         isl_dim_copy(fold->dim), fold->n);
2642         if (!dup)
2643                 return NULL;
2644
2645         for (i = 0; i < fold->n; ++i) {
2646                 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
2647                 if (!dup->qp[i])
2648                         goto error;
2649         }
2650
2651         return dup;
2652 error:
2653         isl_qpolynomial_fold_free(dup);
2654         return NULL;
2655 }
2656
2657 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
2658         __isl_take isl_qpolynomial_fold *fold)
2659 {
2660         if (!fold)
2661                 return NULL;
2662
2663         if (fold->ref == 1)
2664                 return fold;
2665         fold->ref--;
2666         return isl_qpolynomial_fold_dup(fold);
2667 }
2668
2669 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
2670 {
2671         int i;
2672
2673         if (!fold)
2674                 return;
2675         if (--fold->ref > 0)
2676                 return;
2677
2678         for (i = 0; i < fold->n; ++i)
2679                 isl_qpolynomial_free(fold->qp[i]);
2680         isl_dim_free(fold->dim);
2681         free(fold);
2682 }
2683
2684 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
2685 {
2686         if (!fold)
2687                 return -1;
2688
2689         return fold->n == 0;
2690 }
2691
2692 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
2693         __isl_take isl_qpolynomial_fold *fold1,
2694         __isl_take isl_qpolynomial_fold *fold2)
2695 {
2696         int i;
2697         struct isl_qpolynomial_fold *res = NULL;
2698
2699         if (!fold1 || !fold2)
2700                 goto error;
2701
2702         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
2703         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
2704                         goto error);
2705
2706         if (isl_qpolynomial_fold_is_empty(fold1)) {
2707                 isl_qpolynomial_fold_free(fold1);
2708                 return fold2;
2709         }
2710
2711         if (isl_qpolynomial_fold_is_empty(fold2)) {
2712                 isl_qpolynomial_fold_free(fold2);
2713                 return fold1;
2714         }
2715
2716         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
2717                                         fold1->n + fold2->n);
2718         if (!res)
2719                 goto error;
2720
2721         for (i = 0; i < fold1->n; ++i) {
2722                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
2723                 if (!res->qp[res->n])
2724                         goto error;
2725                 res->n++;
2726         }
2727
2728         for (i = 0; i < fold2->n; ++i) {
2729                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
2730                 if (!res->qp[res->n])
2731                         goto error;
2732                 res->n++;
2733         }
2734
2735         isl_qpolynomial_fold_free(fold1);
2736         isl_qpolynomial_fold_free(fold2);
2737
2738         return res;
2739 error:
2740         isl_qpolynomial_fold_free(res);
2741         isl_qpolynomial_fold_free(fold1);
2742         isl_qpolynomial_fold_free(fold2);
2743         return NULL;
2744 }
2745
2746 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2747         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
2748 {
2749         int i;
2750         isl_pw_qpolynomial_fold *pwf;
2751
2752         if (!pwqp)
2753                 return NULL;
2754
2755         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
2756
2757         for (i = 0; i < pwqp->n; ++i)
2758                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
2759                         isl_set_copy(pwqp->p[i].set),
2760                         isl_qpolynomial_fold_alloc(type,
2761                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
2762
2763         isl_pw_qpolynomial_free(pwqp);
2764
2765         return pwf;
2766 }
2767
2768 /* For each parameter or variable that does not appear in qp,
2769  * first eliminate the variable from all constraints and then set it to zero.
2770  */
2771 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
2772         __isl_keep isl_qpolynomial *qp)
2773 {
2774         int *active = NULL;
2775         int i;
2776         int d;
2777         unsigned nparam;
2778         unsigned nvar;
2779
2780         if (!set || !qp)
2781                 goto error;
2782
2783         d = isl_dim_total(set->dim);
2784         active = isl_calloc_array(set->ctx, int, d);
2785         if (set_active(qp, active) < 0)
2786                 goto error;
2787
2788         for (i = 0; i < d; ++i)
2789                 if (!active[i])
2790                         break;
2791
2792         if (i == d) {
2793                 free(active);
2794                 return set;
2795         }
2796
2797         nparam = isl_dim_size(set->dim, isl_dim_param);
2798         nvar = isl_dim_size(set->dim, isl_dim_set);
2799         for (i = 0; i < nparam; ++i) {
2800                 if (active[i])
2801                         continue;
2802                 set = isl_set_eliminate(set, isl_dim_param, i, 1);
2803                 set = isl_set_fix_si(set, isl_dim_param, i, 0);
2804         }
2805         for (i = 0; i < nvar; ++i) {
2806                 if (active[i])
2807                         continue;
2808                 set = isl_set_eliminate(set, isl_dim_set, i, 1);
2809                 set = isl_set_fix_si(set, isl_dim_set, i, 0);
2810         }
2811
2812         free(active);
2813
2814         return set;
2815 error:
2816         free(active);
2817         isl_set_free(set);
2818         return NULL;
2819 }
2820
2821 struct isl_max_data {
2822         isl_qpolynomial *qp;
2823         int first;
2824         isl_qpolynomial *max;
2825 };
2826
2827 static int max_fn(__isl_take isl_point *pnt, void *user)
2828 {
2829         struct isl_max_data *data = (struct isl_max_data *)user;
2830         isl_qpolynomial *val;
2831
2832         val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
2833         if (data->first) {
2834                 data->first = 0;
2835                 data->max = val;
2836         } else {
2837                 data->max = qpolynomial_max(data->max, val);
2838         }
2839
2840         return 0;
2841 }
2842
2843 static __isl_give isl_qpolynomial *guarded_qpolynomial_max(
2844         __isl_take isl_set *set, __isl_take isl_qpolynomial *qp)
2845 {
2846         struct isl_max_data data = { NULL, 1, NULL };
2847
2848         if (!set || !qp)
2849                 goto error;
2850
2851         if (isl_upoly_is_cst(qp->upoly)) {
2852                 isl_set_free(set);
2853                 return qp;
2854         }
2855
2856         set = fix_inactive(set, qp);
2857
2858         data.qp = qp;
2859         if (isl_set_foreach_point(set, max_fn, &data) < 0)
2860                 goto error;
2861
2862         isl_set_free(set);
2863         isl_qpolynomial_free(qp);
2864         return data.max;
2865 error:
2866         isl_set_free(set);
2867         isl_qpolynomial_free(qp);
2868         isl_qpolynomial_free(data.max);
2869         return NULL;
2870 }
2871
2872 /* Compute the maximal value attained by the piecewise quasipolynomial
2873  * on its domain or zero if the domain is empty.
2874  * In the worst case, the domain is scanned completely,
2875  * so the domain is assumed to be bounded.
2876  */
2877 __isl_give isl_qpolynomial *isl_pw_qpolynomial_max(
2878         __isl_take isl_pw_qpolynomial *pwqp)
2879 {
2880         int i;
2881         isl_qpolynomial *max;
2882
2883         if (!pwqp)
2884                 return NULL;
2885
2886         if (pwqp->n == 0) {
2887                 isl_dim *dim = isl_dim_copy(pwqp->dim);
2888                 isl_pw_qpolynomial_free(pwqp);
2889                 return isl_qpolynomial_zero(dim);
2890         }
2891
2892         max = guarded_qpolynomial_max(isl_set_copy(pwqp->p[0].set),
2893                                         isl_qpolynomial_copy(pwqp->p[0].qp));
2894         for (i = 1; i < pwqp->n; ++i) {
2895                 isl_qpolynomial *max_i;
2896                 max_i = guarded_qpolynomial_max(isl_set_copy(pwqp->p[i].set),
2897                                             isl_qpolynomial_copy(pwqp->p[i].qp));
2898                 max = qpolynomial_max(max, max_i);
2899         }
2900
2901         isl_pw_qpolynomial_free(pwqp);
2902         return max;
2903 }