isl_polynomial.c

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