isl_polynomial.c

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