isl_fold.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 <isl_union_map_private.h>
  12 #include <isl_polynomial_private.h>
  13 #include <isl_point_private.h>
  14 #include <isl_dim_private.h>
  15 #include <isl_map_private.h>
  16 #include <isl_lp.h>
  17 #include <isl_seq.h>
  18
  19 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
  20         enum isl_fold type, __isl_take isl_dim *dim, int n)
  21 {
  22         isl_qpolynomial_fold *fold;
  23
  24         if (!dim)
  25                 goto error;
  26
  27         isl_assert(dim->ctx, n >= 0, goto error);
  28         fold = isl_calloc(dim->ctx, struct isl_qpolynomial_fold,
  29                         sizeof(struct isl_qpolynomial_fold) +
  30                         (n - 1) * sizeof(struct isl_qpolynomial *));
  31         if (!fold)
  32                 goto error;
  33
  34         fold->ref = 1;
  35         fold->size = n;
  36         fold->n = 0;
  37         fold->type = type;
  38         fold->dim = dim;
  39
  40         return fold;
  41 error:
  42         isl_dim_free(dim);
  43         return NULL;
  44 }
  45
  46 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_reset_dim(
  47         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_dim *dim)
  48 {
  49         int i;
  50
  51         fold = isl_qpolynomial_fold_cow(fold);
  52         if (!fold || !dim)
  53                 goto error;
  54
  55         for (i = 0; i < fold->n; ++i) {
  56                 fold->qp[i] = isl_qpolynomial_reset_dim(fold->qp[i],
  57                                                         isl_dim_copy(dim));
  58                 if (!fold->qp[i])
  59                         goto error;
  60         }
  61
  62         isl_dim_free(fold->dim);
  63         fold->dim = dim;
  64
  65         return fold;
  66 error:
  67         isl_qpolynomial_fold_free(fold);
  68         isl_dim_free(dim);
  69         return NULL;
  70 }
  71
  72 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
  73         enum isl_dim_type type, unsigned first, unsigned n)
  74 {
  75         int i;
  76
  77         if (!fold)
  78                 return -1;
  79         if (fold->n == 0 || n == 0)
  80                 return 0;
  81
  82         for (i = 0; i < fold->n; ++i) {
  83                 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
  84                                                             type, first, n);
  85                 if (involves < 0 || involves)
  86                         return involves;
  87         }
  88         return 0;
  89 }
  90
  91 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
  92         __isl_take isl_qpolynomial_fold *fold,
  93         enum isl_dim_type type, unsigned first, unsigned n)
  94 {
  95         int i;
  96
  97         if (!fold)
  98                 return NULL;
  99         if (n == 0)
 100                 return fold;
 101
 102         fold = isl_qpolynomial_fold_cow(fold);
 103         if (!fold)
 104                 return NULL;
 105         fold->dim = isl_dim_drop(fold->dim, type, first, n);
 106         if (!fold->dim)
 107                 goto error;
 108
 109         for (i = 0; i < fold->n; ++i) {
 110                 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
 111                                                             type, first, n);
 112                 if (!fold->qp[i])
 113                         goto error;
 114         }
 115
 116         return fold;
 117 error:
 118         isl_qpolynomial_fold_free(fold);
 119         return NULL;
 120 }
 121
 122 static int isl_qpolynomial_cst_sign(__isl_keep isl_qpolynomial *qp)
 123 {
 124         struct isl_upoly_cst *cst;
 125
 126         cst = isl_upoly_as_cst(qp->upoly);
 127         if (!cst)
 128                 return 0;
 129
 130         return isl_int_sgn(cst->n) < 0 ? -1 : 1;
 131 }
 132
 133 static int isl_qpolynomial_aff_sign(__isl_keep isl_set *set,
 134         __isl_keep isl_qpolynomial *qp)
 135 {
 136         enum isl_lp_result res;
 137         isl_vec *aff;
 138         isl_int opt;
 139         int sgn = 0;
 140
 141         aff = isl_qpolynomial_extract_affine(qp);
 142         if (!aff)
 143                 return 0;
 144
 145         isl_int_init(opt);
 146
 147         res = isl_set_solve_lp(set, 0, aff->el + 1, aff->el[0],
 148                                 &opt, NULL, NULL);
 149         if (res == isl_lp_error)
 150                 goto done;
 151         if (res == isl_lp_empty ||
 152             (res == isl_lp_ok && !isl_int_is_neg(opt))) {
 153                 sgn = 1;
 154                 goto done;
 155         }
 156
 157         res = isl_set_solve_lp(set, 1, aff->el + 1, aff->el[0],
 158                                 &opt, NULL, NULL);
 159         if (res == isl_lp_ok && !isl_int_is_pos(opt))
 160                 sgn = -1;
 161
 162 done:
 163         isl_int_clear(opt);
 164         isl_vec_free(aff);
 165         return sgn;
 166 }
 167
 168 /* Determine, if possible, the sign of the quasipolynomial "qp" on
 169  * the domain "set".
 170  *
 171  * If qp is a constant, then the problem is trivial.
 172  * If qp is linear, then we check if the minimum of the corresponding
 173  * affine constraint is non-negative or if the maximum is non-positive.
 174  *
 175  * Otherwise, we check if the outermost variable "v" has a lower bound "l"
 176  * in "set".  If so, we write qp(v,v') as
 177  *
 178  *      q(v,v') * (v - l) + r(v')
 179  *
 180  * if q(v,v') and r(v') have the same known sign, then the original
 181  * quasipolynomial has the same sign as well.
 182  *
 183  * Return
 184  *      -1 if qp <= 0
 185  *       1 if qp >= 0
 186  *       0 if unknown
 187  */
 188 static int isl_qpolynomial_sign(__isl_keep isl_set *set,
 189         __isl_keep isl_qpolynomial *qp)
 190 {
 191         int d;
 192         int i;
 193         int is;
 194         struct isl_upoly_rec *rec;
 195         isl_vec *v;
 196         isl_int l;
 197         enum isl_lp_result res;
 198         int sgn = 0;
 199
 200         is = isl_qpolynomial_is_cst(qp, NULL, NULL);
 201         if (is < 0)
 202                 return 0;
 203         if (is)
 204                 return isl_qpolynomial_cst_sign(qp);
 205
 206         is = isl_qpolynomial_is_affine(qp);
 207         if (is < 0)
 208                 return 0;
 209         if (is)
 210                 return isl_qpolynomial_aff_sign(set, qp);
 211
 212         if (qp->div->n_row > 0)
 213                 return 0;
 214
 215         rec = isl_upoly_as_rec(qp->upoly);
 216         if (!rec)
 217                 return 0;
 218
 219         d = isl_dim_total(qp->dim);
 220         v = isl_vec_alloc(set->ctx, 2 + d);
 221         if (!v)
 222                 return 0;
 223
 224         isl_seq_clr(v->el + 1, 1 + d);
 225         isl_int_set_si(v->el[0], 1);
 226         isl_int_set_si(v->el[2 + qp->upoly->var], 1);
 227
 228         isl_int_init(l);
 229
 230         res = isl_set_solve_lp(set, 0, v->el + 1, v->el[0], &l, NULL, NULL);
 231         if (res == isl_lp_ok) {
 232                 isl_qpolynomial *min;
 233                 isl_qpolynomial *base;
 234                 isl_qpolynomial *r, *q;
 235                 isl_qpolynomial *t;
 236
 237                 min = isl_qpolynomial_cst(isl_dim_copy(qp->dim), l);
 238                 base = isl_qpolynomial_pow(isl_dim_copy(qp->dim),
 239                                                 qp->upoly->var, 1);
 240
 241                 r = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
 242                                           isl_upoly_copy(rec->p[rec->n - 1]));
 243                 q = isl_qpolynomial_copy(r);
 244
 245                 for (i = rec->n - 2; i >= 0; --i) {
 246                         r = isl_qpolynomial_mul(r, isl_qpolynomial_copy(min));
 247                         t = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
 248                                                   isl_upoly_copy(rec->p[i]));
 249                         r = isl_qpolynomial_add(r, t);
 250                         if (i == 0)
 251                                 break;
 252                         q = isl_qpolynomial_mul(q, isl_qpolynomial_copy(base));
 253                         q = isl_qpolynomial_add(q, isl_qpolynomial_copy(r));
 254                 }
 255
 256                 if (isl_qpolynomial_is_zero(q))
 257                         sgn = isl_qpolynomial_sign(set, r);
 258                 else if (isl_qpolynomial_is_zero(r))
 259                         sgn = isl_qpolynomial_sign(set, q);
 260                 else {
 261                         int sgn_q, sgn_r;
 262                         sgn_r = isl_qpolynomial_sign(set, r);
 263                         sgn_q = isl_qpolynomial_sign(set, q);
 264                         if (sgn_r == sgn_q)
 265                                 sgn = sgn_r;
 266                 }
 267
 268                 isl_qpolynomial_free(min);
 269                 isl_qpolynomial_free(base);
 270                 isl_qpolynomial_free(q);
 271                 isl_qpolynomial_free(r);
 272         }
 273
 274         isl_int_clear(l);
 275
 276         isl_vec_free(v);
 277
 278         return sgn;
 279 }
 280
 281 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold_on_domain(
 282         __isl_keep isl_set *set,
 283         __isl_take isl_qpolynomial_fold *fold1,
 284         __isl_take isl_qpolynomial_fold *fold2)
 285 {
 286         int i, j;
 287         int n1;
 288         struct isl_qpolynomial_fold *res = NULL;
 289         int better;
 290
 291         if (!fold1 || !fold2)
 292                 goto error;
 293
 294         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
 295         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
 296                         goto error);
 297
 298         better = fold1->type == isl_fold_max ? -1 : 1;
 299
 300         if (isl_qpolynomial_fold_is_empty(fold1)) {
 301                 isl_qpolynomial_fold_free(fold1);
 302                 return fold2;
 303         }
 304
 305         if (isl_qpolynomial_fold_is_empty(fold2)) {
 306                 isl_qpolynomial_fold_free(fold2);
 307                 return fold1;
 308         }
 309
 310         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
 311                                         fold1->n + fold2->n);
 312         if (!res)
 313                 goto error;
 314
 315         for (i = 0; i < fold1->n; ++i) {
 316                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
 317                 if (!res->qp[res->n])
 318                         goto error;
 319                 res->n++;
 320         }
 321         n1 = res->n;
 322
 323         for (i = 0; i < fold2->n; ++i) {
 324                 for (j = n1 - 1; j >= 0; --j) {
 325                         isl_qpolynomial *d;
 326                         int sgn;
 327                         d = isl_qpolynomial_sub(
 328                                 isl_qpolynomial_copy(res->qp[j]),
 329                                 isl_qpolynomial_copy(fold2->qp[i]));
 330                         sgn = isl_qpolynomial_sign(set, d);
 331                         isl_qpolynomial_free(d);
 332                         if (sgn == 0)
 333                                 continue;
 334                         if (sgn != better)
 335                                 break;
 336                         isl_qpolynomial_free(res->qp[j]);
 337                         if (j != n1 - 1)
 338                                 res->qp[j] = res->qp[n1 - 1];
 339                         n1--;
 340                         if (n1 != res->n - 1)
 341                                 res->qp[n1] = res->qp[res->n - 1];
 342                         res->n--;
 343                 }
 344                 if (j >= 0)
 345                         continue;
 346                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
 347                 if (!res->qp[res->n])
 348                         goto error;
 349                 res->n++;
 350         }
 351
 352         isl_qpolynomial_fold_free(fold1);
 353         isl_qpolynomial_fold_free(fold2);
 354
 355         return res;
 356 error:
 357         isl_qpolynomial_fold_free(res);
 358         isl_qpolynomial_fold_free(fold1);
 359         isl_qpolynomial_fold_free(fold2);
 360         return NULL;
 361 }
 362
 363 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_add_qpolynomial(
 364         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_qpolynomial *qp)
 365 {
 366         int i;
 367
 368         if (!fold || !qp)
 369                 goto error;
 370
 371         if (isl_qpolynomial_is_zero(qp)) {
 372                 isl_qpolynomial_free(qp);
 373                 return fold;
 374         }
 375
 376         fold = isl_qpolynomial_fold_cow(fold);
 377         if (!fold)
 378                 goto error;
 379
 380         for (i = 0; i < fold->n; ++i) {
 381                 fold->qp[i] = isl_qpolynomial_add(fold->qp[i],
 382                                                 isl_qpolynomial_copy(qp));
 383                 if (!fold->qp[i])
 384                         goto error;
 385         }
 386
 387         isl_qpolynomial_free(qp);
 388         return fold;
 389 error:
 390         isl_qpolynomial_fold_free(fold);
 391         isl_qpolynomial_free(qp);
 392         return NULL;
 393 }
 394
 395 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_add_on_domain(
 396         __isl_keep isl_set *dom,
 397         __isl_take isl_qpolynomial_fold *fold1,
 398         __isl_take isl_qpolynomial_fold *fold2)
 399 {
 400         int i;
 401         isl_qpolynomial_fold *res = NULL;
 402
 403         if (!fold1 || !fold2)
 404                 goto error;
 405
 406         if (isl_qpolynomial_fold_is_empty(fold1)) {
 407                 isl_qpolynomial_fold_free(fold1);
 408                 return fold2;
 409         }
 410
 411         if (isl_qpolynomial_fold_is_empty(fold2)) {
 412                 isl_qpolynomial_fold_free(fold2);
 413                 return fold1;
 414         }
 415
 416         if (fold1->n == 1 && fold2->n != 1)
 417                 return isl_qpolynomial_fold_add_on_domain(dom, fold2, fold1);
 418
 419         if (fold2->n == 1) {
 420                 res = isl_qpolynomial_fold_add_qpolynomial(fold1,
 421                                         isl_qpolynomial_copy(fold2->qp[0]));
 422                 isl_qpolynomial_fold_free(fold2);
 423                 return res;
 424         }
 425
 426         res = isl_qpolynomial_fold_add_qpolynomial(
 427                                 isl_qpolynomial_fold_copy(fold1),
 428                                 isl_qpolynomial_copy(fold2->qp[0]));
 429
 430         for (i = 1; i < fold2->n; ++i) {
 431                 isl_qpolynomial_fold *res_i;
 432                 res_i = isl_qpolynomial_fold_add_qpolynomial(
 433                                         isl_qpolynomial_fold_copy(fold1),
 434                                         isl_qpolynomial_copy(fold2->qp[i]));
 435                 res = isl_qpolynomial_fold_fold_on_domain(dom, res, res_i);
 436         }
 437
 438         isl_qpolynomial_fold_free(fold1);
 439         isl_qpolynomial_fold_free(fold2);
 440         return res;
 441 error:
 442         isl_qpolynomial_fold_free(res);
 443         isl_qpolynomial_fold_free(fold1);
 444         isl_qpolynomial_fold_free(fold2);
 445         return NULL;
 446 }
 447
 448 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_substitute_equalities(
 449         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_basic_set *eq)
 450 {
 451         int i;
 452
 453         if (!fold || !eq)
 454                 goto error;
 455
 456         fold = isl_qpolynomial_fold_cow(fold);
 457         if (!fold)
 458                 return NULL;
 459
 460         for (i = 0; i < fold->n; ++i) {
 461                 fold->qp[i] = isl_qpolynomial_substitute_equalities(fold->qp[i],
 462                                                         isl_basic_set_copy(eq));
 463                 if (!fold->qp[i])
 464                         goto error;
 465         }
 466
 467         isl_basic_set_free(eq);
 468         return fold;
 469 error:
 470         isl_basic_set_free(eq);
 471         isl_qpolynomial_fold_free(fold);
 472         return NULL;
 473 }
 474
 475 #define HAS_TYPE
 476
 477 #undef PW
 478 #define PW isl_pw_qpolynomial_fold
 479 #undef EL
 480 #define EL isl_qpolynomial_fold
 481 #undef IS_ZERO
 482 #define IS_ZERO is_empty
 483 #undef FIELD
 484 #define FIELD fold
 485
 486 #include <isl_pw_templ.c>
 487
 488 #undef UNION
 489 #define UNION isl_union_pw_qpolynomial_fold
 490 #undef PART
 491 #define PART isl_pw_qpolynomial_fold
 492 #undef PARTS
 493 #define PARTS pw_qpolynomial_fold
 494
 495 #include <isl_union_templ.c>
 496
 497 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
 498         __isl_take isl_dim *dim)
 499 {
 500         return qpolynomial_fold_alloc(type, dim, 0);
 501 }
 502
 503 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
 504         enum isl_fold type, __isl_take isl_qpolynomial *qp)
 505 {
 506         isl_qpolynomial_fold *fold;
 507
 508         if (!qp)
 509                 return NULL;
 510
 511         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
 512         if (!fold)
 513                 goto error;
 514
 515         fold->qp[0] = qp;
 516         fold->n++;
 517
 518         return fold;
 519 error:
 520         isl_qpolynomial_fold_free(fold);
 521         isl_qpolynomial_free(qp);
 522         return NULL;
 523 }
 524
 525 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
 526         __isl_keep isl_qpolynomial_fold *fold)
 527 {
 528         if (!fold)
 529                 return NULL;
 530
 531         fold->ref++;
 532         return fold;
 533 }
 534
 535 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
 536         __isl_keep isl_qpolynomial_fold *fold)
 537 {
 538         int i;
 539         isl_qpolynomial_fold *dup;
 540
 541         if (!fold)
 542                 return NULL;
 543         dup = qpolynomial_fold_alloc(fold->type,
 544                                         isl_dim_copy(fold->dim), fold->n);
 545         if (!dup)
 546                 return NULL;
 547
 548         dup->n = fold->n;
 549         for (i = 0; i < fold->n; ++i) {
 550                 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
 551                 if (!dup->qp[i])
 552                         goto error;
 553         }
 554
 555         return dup;
 556 error:
 557         isl_qpolynomial_fold_free(dup);
 558         return NULL;
 559 }
 560
 561 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
 562         __isl_take isl_qpolynomial_fold *fold)
 563 {
 564         if (!fold)
 565                 return NULL;
 566
 567         if (fold->ref == 1)
 568                 return fold;
 569         fold->ref--;
 570         return isl_qpolynomial_fold_dup(fold);
 571 }
 572
 573 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
 574 {
 575         int i;
 576
 577         if (!fold)
 578                 return;
 579         if (--fold->ref > 0)
 580                 return;
 581
 582         for (i = 0; i < fold->n; ++i)
 583                 isl_qpolynomial_free(fold->qp[i]);
 584         isl_dim_free(fold->dim);
 585         free(fold);
 586 }
 587
 588 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
 589 {
 590         if (!fold)
 591                 return -1;
 592
 593         return fold->n == 0;
 594 }
 595
 596 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
 597         __isl_take isl_qpolynomial_fold *fold1,
 598         __isl_take isl_qpolynomial_fold *fold2)
 599 {
 600         int i;
 601         struct isl_qpolynomial_fold *res = NULL;
 602
 603         if (!fold1 || !fold2)
 604                 goto error;
 605
 606         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
 607         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
 608                         goto error);
 609
 610         if (isl_qpolynomial_fold_is_empty(fold1)) {
 611                 isl_qpolynomial_fold_free(fold1);
 612                 return fold2;
 613         }
 614
 615         if (isl_qpolynomial_fold_is_empty(fold2)) {
 616                 isl_qpolynomial_fold_free(fold2);
 617                 return fold1;
 618         }
 619
 620         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
 621                                         fold1->n + fold2->n);
 622         if (!res)
 623                 goto error;
 624
 625         for (i = 0; i < fold1->n; ++i) {
 626                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
 627                 if (!res->qp[res->n])
 628                         goto error;
 629                 res->n++;
 630         }
 631
 632         for (i = 0; i < fold2->n; ++i) {
 633                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
 634                 if (!res->qp[res->n])
 635                         goto error;
 636                 res->n++;
 637         }
 638
 639         isl_qpolynomial_fold_free(fold1);
 640         isl_qpolynomial_fold_free(fold2);
 641
 642         return res;
 643 error:
 644         isl_qpolynomial_fold_free(res);
 645         isl_qpolynomial_fold_free(fold1);
 646         isl_qpolynomial_fold_free(fold2);
 647         return NULL;
 648 }
 649
 650 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_fold(
 651         __isl_take isl_pw_qpolynomial_fold *pw1,
 652         __isl_take isl_pw_qpolynomial_fold *pw2)
 653 {
 654         int i, j, n;
 655         struct isl_pw_qpolynomial_fold *res;
 656         isl_set *set;
 657
 658         if (!pw1 || !pw2)
 659                 goto error;
 660
 661         isl_assert(pw1->dim->ctx, isl_dim_equal(pw1->dim, pw2->dim), goto error);
 662
 663         if (isl_pw_qpolynomial_fold_is_zero(pw1)) {
 664                 isl_pw_qpolynomial_fold_free(pw1);
 665                 return pw2;
 666         }
 667
 668         if (isl_pw_qpolynomial_fold_is_zero(pw2)) {
 669                 isl_pw_qpolynomial_fold_free(pw2);
 670                 return pw1;
 671         }
 672
 673         if (pw1->type != pw2->type)
 674                 isl_die(set->ctx, isl_error_invalid, "fold types don't match",
 675                         goto error);
 676
 677         n = (pw1->n + 1) * (pw2->n + 1);
 678         res = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pw1->dim),
 679                                                 pw1->type, n);
 680
 681         for (i = 0; i < pw1->n; ++i) {
 682                 set = isl_set_copy(pw1->p[i].set);
 683                 for (j = 0; j < pw2->n; ++j) {
 684                         struct isl_set *common;
 685                         isl_qpolynomial_fold *sum;
 686                         set = isl_set_subtract(set,
 687                                         isl_set_copy(pw2->p[j].set));
 688                         common = isl_set_intersect(isl_set_copy(pw1->p[i].set),
 689                                                 isl_set_copy(pw2->p[j].set));
 690                         if (isl_set_fast_is_empty(common)) {
 691                                 isl_set_free(common);
 692                                 continue;
 693                         }
 694
 695                         sum = isl_qpolynomial_fold_fold_on_domain(common,
 696                                isl_qpolynomial_fold_copy(pw1->p[i].fold),
 697                                isl_qpolynomial_fold_copy(pw2->p[j].fold));
 698
 699                         res = isl_pw_qpolynomial_fold_add_piece(res, common, sum);
 700                 }
 701                 res = isl_pw_qpolynomial_fold_add_piece(res, set,
 702                         isl_qpolynomial_fold_copy(pw1->p[i].fold));
 703         }
 704
 705         for (j = 0; j < pw2->n; ++j) {
 706                 set = isl_set_copy(pw2->p[j].set);
 707                 for (i = 0; i < pw1->n; ++i)
 708                         set = isl_set_subtract(set, isl_set_copy(pw1->p[i].set));
 709                 res = isl_pw_qpolynomial_fold_add_piece(res, set,
 710                                     isl_qpolynomial_fold_copy(pw2->p[j].fold));
 711         }
 712
 713         isl_pw_qpolynomial_fold_free(pw1);
 714         isl_pw_qpolynomial_fold_free(pw2);
 715
 716         return res;
 717 error:
 718         isl_pw_qpolynomial_fold_free(pw1);
 719         isl_pw_qpolynomial_fold_free(pw2);
 720         return NULL;
 721 }
 722
 723 __isl_give isl_union_pw_qpolynomial_fold *isl_union_pw_qpolynomial_fold_fold_pw_qpolynomial_fold(
 724         __isl_take isl_union_pw_qpolynomial_fold *u,
 725         __isl_take isl_pw_qpolynomial_fold *part)
 726 {
 727         uint32_t hash;
 728         struct isl_hash_table_entry *entry;
 729
 730         u = isl_union_pw_qpolynomial_fold_cow(u);
 731
 732         if (!part || !u)
 733                 goto error;
 734
 735         isl_assert(u->dim->ctx, isl_dim_match(part->dim, isl_dim_param, u->dim,
 736                                               isl_dim_param), goto error);
 737
 738         hash = isl_dim_get_hash(part->dim);
 739         entry = isl_hash_table_find(u->dim->ctx, &u->table, hash,
 740                                     &has_dim, part->dim, 1);
 741         if (!entry)
 742                 goto error;
 743
 744         if (!entry->data)
 745                 entry->data = part;
 746         else {
 747                 entry->data = isl_pw_qpolynomial_fold_fold(entry->data,
 748                                             isl_pw_qpolynomial_fold_copy(part));
 749                 if (!entry->data)
 750                         goto error;
 751                 isl_pw_qpolynomial_fold_free(part);
 752         }
 753
 754         return u;
 755 error:
 756         isl_pw_qpolynomial_fold_free(part);
 757         isl_union_pw_qpolynomial_fold_free(u);
 758         return NULL;
 759 }
 760
 761 static int fold_part(__isl_take isl_pw_qpolynomial_fold *part, void *user)
 762 {
 763         isl_union_pw_qpolynomial_fold **u;
 764         u = (isl_union_pw_qpolynomial_fold **)user;
 765
 766         *u = isl_union_pw_qpolynomial_fold_fold_pw_qpolynomial_fold(*u, part);
 767
 768         return 0;
 769 }
 770
 771 __isl_give isl_union_pw_qpolynomial_fold *isl_union_pw_qpolynomial_fold_fold(
 772         __isl_take isl_union_pw_qpolynomial_fold *u1,
 773         __isl_take isl_union_pw_qpolynomial_fold *u2)
 774 {
 775         u1 = isl_union_pw_qpolynomial_fold_cow(u1);
 776
 777         if (!u1 || !u2)
 778                 goto error;
 779
 780         if (isl_union_pw_qpolynomial_fold_foreach_pw_qpolynomial_fold(u2,
 781                                                         &fold_part, &u1) < 0)
 782                 goto error;
 783
 784         isl_union_pw_qpolynomial_fold_free(u2);
 785
 786         return u1;
 787 error:
 788         isl_union_pw_qpolynomial_fold_free(u1);
 789         isl_union_pw_qpolynomial_fold_free(u2);
 790         return NULL;
 791 }
 792
 793 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
 794         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
 795 {
 796         int i;
 797         isl_pw_qpolynomial_fold *pwf;
 798
 799         if (!pwqp)
 800                 return NULL;
 801
 802         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), type,
 803                                                 pwqp->n);
 804
 805         for (i = 0; i < pwqp->n; ++i)
 806                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
 807                         isl_set_copy(pwqp->p[i].set),
 808                         isl_qpolynomial_fold_alloc(type,
 809                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
 810
 811         isl_pw_qpolynomial_free(pwqp);
 812
 813         return pwf;
 814 }
 815
 816 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
 817         __isl_keep isl_qpolynomial_fold *fold2)
 818 {
 819         int i;
 820
 821         if (!fold1 || !fold2)
 822                 return -1;
 823
 824         if (fold1->n != fold2->n)
 825                 return 0;
 826
 827         /* We probably want to sort the qps first... */
 828         for (i = 0; i < fold1->n; ++i) {
 829                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
 830                 if (eq < 0 || !eq)
 831                         return eq;
 832         }
 833
 834         return 1;
 835 }
 836
 837 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
 838         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
 839 {
 840         isl_qpolynomial *qp;
 841
 842         if (!fold || !pnt)
 843                 goto error;
 844         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
 845         isl_assert(pnt->dim->ctx,
 846                 fold->type == isl_fold_max || fold->type == isl_fold_min,
 847                 goto error);
 848
 849         if (fold->n == 0)
 850                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
 851         else {
 852                 int i;
 853                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
 854                                                 isl_point_copy(pnt));
 855                 for (i = 1; i < fold->n; ++i) {
 856                         isl_qpolynomial *qp_i;
 857                         qp_i = isl_qpolynomial_eval(
 858                                             isl_qpolynomial_copy(fold->qp[i]),
 859                                             isl_point_copy(pnt));
 860                         if (fold->type == isl_fold_max)
 861                                 qp = isl_qpolynomial_max_cst(qp, qp_i);
 862                         else
 863                                 qp = isl_qpolynomial_min_cst(qp, qp_i);
 864                 }
 865         }
 866         isl_qpolynomial_fold_free(fold);
 867         isl_point_free(pnt);
 868
 869         return qp;
 870 error:
 871         isl_qpolynomial_fold_free(fold);
 872         isl_point_free(pnt);
 873         return NULL;
 874 }
 875
 876 size_t isl_pw_qpolynomial_fold_size(__isl_keep isl_pw_qpolynomial_fold *pwf)
 877 {
 878         int i;
 879         size_t n = 0;
 880
 881         for (i = 0; i < pwf->n; ++i)
 882                 n += pwf->p[i].fold->n;
 883
 884         return n;
 885 }
 886
 887 __isl_give isl_qpolynomial *isl_qpolynomial_fold_opt_on_domain(
 888         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *set, int max)
 889 {
 890         int i;
 891         isl_qpolynomial *opt;
 892
 893         if (!set || !fold)
 894                 goto error;
 895
 896         if (fold->n == 0) {
 897                 isl_dim *dim = isl_dim_copy(fold->dim);
 898                 isl_set_free(set);
 899                 isl_qpolynomial_fold_free(fold);
 900                 return isl_qpolynomial_zero(dim);
 901         }
 902
 903         opt = isl_qpolynomial_opt_on_domain(isl_qpolynomial_copy(fold->qp[0]),
 904                                                 isl_set_copy(set), max);
 905         for (i = 1; i < fold->n; ++i) {
 906                 isl_qpolynomial *opt_i;
 907                 opt_i = isl_qpolynomial_opt_on_domain(
 908                                 isl_qpolynomial_copy(fold->qp[i]),
 909                                 isl_set_copy(set), max);
 910                 if (max)
 911                         opt = isl_qpolynomial_max_cst(opt, opt_i);
 912                 else
 913                         opt = isl_qpolynomial_min_cst(opt, opt_i);
 914         }
 915
 916         isl_set_free(set);
 917         isl_qpolynomial_fold_free(fold);
 918
 919         return opt;
 920 error:
 921         isl_set_free(set);
 922         isl_qpolynomial_fold_free(fold);
 923         return NULL;
 924 }
 925
 926 /* Check whether for each quasi-polynomial in "fold2" there is
 927  * a quasi-polynomial in "fold1" that dominates it on "set".
 928  */
 929 static int qpolynomial_fold_covers_on_domain(__isl_keep isl_set *set,
 930         __isl_keep isl_qpolynomial_fold *fold1,
 931         __isl_keep isl_qpolynomial_fold *fold2)
 932 {
 933         int i, j;
 934         int covers;
 935
 936         if (!set || !fold1 || !fold2)
 937                 return -1;
 938
 939         covers = fold1->type == isl_fold_max ? 1 : -1;
 940
 941         for (i = 0; i < fold2->n; ++i) {
 942                 for (j = 0; j < fold1->n; ++j) {
 943                         isl_qpolynomial *d;
 944                         int sgn;
 945
 946                         d = isl_qpolynomial_sub(
 947                                 isl_qpolynomial_copy(fold1->qp[j]),
 948                                 isl_qpolynomial_copy(fold2->qp[i]));
 949                         sgn = isl_qpolynomial_sign(set, d);
 950                         isl_qpolynomial_free(d);
 951                         if (sgn == covers)
 952                                 break;
 953                 }
 954                 if (j >= fold1->n)
 955                         return 0;
 956         }
 957
 958         return 1;
 959 }
 960
 961 /* Check whether "pwf1" dominated "pwf2", i.e., the domain of "pwf1" contains
 962  * that of "pwf2" and on each cell, the corresponding fold from pwf1 dominates
 963  * that of pwf2.
 964  */
 965 int isl_pw_qpolynomial_fold_covers(__isl_keep isl_pw_qpolynomial_fold *pwf1,
 966         __isl_keep isl_pw_qpolynomial_fold *pwf2)
 967 {
 968         int i, j;
 969         isl_set *dom1, *dom2;
 970         int is_subset;
 971
 972         if (!pwf1 || !pwf2)
 973                 return -1;
 974
 975         if (pwf2->n == 0)
 976                 return 1;
 977         if (pwf1->n == 0)
 978                 return 0;
 979
 980         dom1 = isl_pw_qpolynomial_fold_domain(isl_pw_qpolynomial_fold_copy(pwf1));
 981         dom2 = isl_pw_qpolynomial_fold_domain(isl_pw_qpolynomial_fold_copy(pwf2));
 982         is_subset = isl_set_is_subset(dom2, dom1);
 983         isl_set_free(dom1);
 984         isl_set_free(dom2);
 985
 986         if (is_subset < 0 || !is_subset)
 987                 return is_subset;
 988
 989         for (i = 0; i < pwf2->n; ++i) {
 990                 for (j = 0; j < pwf1->n; ++j) {
 991                         int is_empty;
 992                         isl_set *common;
 993                         int covers;
 994
 995                         common = isl_set_intersect(isl_set_copy(pwf1->p[j].set),
 996                                                    isl_set_copy(pwf2->p[i].set));
 997                         is_empty = isl_set_is_empty(common);
 998                         if (is_empty < 0 || is_empty) {
 999                                 isl_set_free(common);
1000                                 if (is_empty < 0)
1001                                         return -1;
1002                                 continue;
1003                         }
1004                         covers = qpolynomial_fold_covers_on_domain(common,
1005                                         pwf1->p[j].fold, pwf2->p[i].fold);
1006                         isl_set_free(common);
1007                         if (covers < 0 || !covers)
1008                                 return covers;
1009                 }
1010         }
1011
1012         return 1;
1013 }
1014
1015 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_morph(
1016         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_morph *morph)
1017 {
1018         int i;
1019         isl_ctx *ctx;
1020
1021         if (!fold || !morph)
1022                 goto error;
1023
1024         ctx = fold->dim->ctx;
1025         isl_assert(ctx, isl_dim_equal(fold->dim, morph->dom->dim), goto error);
1026
1027         fold = isl_qpolynomial_fold_cow(fold);
1028         if (!fold)
1029                 goto error;
1030
1031         isl_dim_free(fold->dim);
1032         fold->dim = isl_dim_copy(morph->ran->dim);
1033         if (!fold->dim)
1034                 goto error;
1035
1036         for (i = 0; i < fold->n; ++i) {
1037                 fold->qp[i] = isl_qpolynomial_morph(fold->qp[i],
1038                                                 isl_morph_copy(morph));
1039                 if (!fold->qp[i])
1040                         goto error;
1041         }
1042
1043         isl_morph_free(morph);
1044
1045         return fold;
1046 error:
1047         isl_qpolynomial_fold_free(fold);
1048         isl_morph_free(morph);
1049         return NULL;
1050 }
1051
1052 enum isl_fold isl_qpolynomial_fold_get_type(__isl_keep isl_qpolynomial_fold *fold)
1053 {
1054         if (!fold)
1055                 return isl_fold_list;
1056         return fold->type;
1057 }
1058
1059 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_lift(
1060         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_dim *dim)
1061 {
1062         int i;
1063         isl_ctx *ctx;
1064
1065         if (!fold || !dim)
1066                 goto error;
1067
1068         if (isl_dim_equal(fold->dim, dim)) {
1069                 isl_dim_free(dim);
1070                 return fold;
1071         }
1072
1073         fold = isl_qpolynomial_fold_cow(fold);
1074         if (!fold)
1075                 goto error;
1076
1077         isl_dim_free(fold->dim);
1078         fold->dim = isl_dim_copy(dim);
1079         if (!fold->dim)
1080                 goto error;
1081
1082         for (i = 0; i < fold->n; ++i) {
1083                 fold->qp[i] = isl_qpolynomial_lift(fold->qp[i],
1084                                                 isl_dim_copy(dim));
1085                 if (!fold->qp[i])
1086                         goto error;
1087         }
1088
1089         isl_dim_free(dim);
1090
1091         return fold;
1092 error:
1093         isl_qpolynomial_fold_free(fold);
1094         isl_dim_free(dim);
1095         return NULL;
1096 }
1097
1098 int isl_qpolynomial_fold_foreach_qpolynomial(
1099         __isl_keep isl_qpolynomial_fold *fold,
1100         int (*fn)(__isl_take isl_qpolynomial *qp, void *user), void *user)
1101 {
1102         int i;
1103
1104         if (!fold)
1105                 return -1;
1106
1107         for (i = 0; i < fold->n; ++i)
1108                 if (fn(isl_qpolynomial_copy(fold->qp[i]), user) < 0)
1109                         return -1;
1110
1111         return 0;
1112 }
1113
1114 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_move_dims(
1115         __isl_take isl_qpolynomial_fold *fold,
1116         enum isl_dim_type dst_type, unsigned dst_pos,
1117         enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1118 {
1119         int i;
1120
1121         if (n == 0)
1122                 return fold;
1123
1124         fold = isl_qpolynomial_fold_cow(fold);
1125         if (!fold)
1126                 return NULL;
1127
1128         fold->dim = isl_dim_move(fold->dim, dst_type, dst_pos,
1129                                                 src_type, src_pos, n);
1130         if (!fold->dim)
1131                 goto error;
1132
1133         for (i = 0; i < fold->n; ++i) {
1134                 fold->qp[i] = isl_qpolynomial_move_dims(fold->qp[i],
1135                                 dst_type, dst_pos, src_type, src_pos, n);
1136                 if (!fold->qp[i])
1137                         goto error;
1138         }
1139
1140         return fold;
1141 error:
1142         isl_qpolynomial_fold_free(fold);
1143         return NULL;
1144 }
1145
1146 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
1147  * in fold->qp[k] by subs[i].
1148  */
1149 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_substitute(
1150         __isl_take isl_qpolynomial_fold *fold,
1151         enum isl_dim_type type, unsigned first, unsigned n,
1152         __isl_keep isl_qpolynomial **subs)
1153 {
1154         int i;
1155
1156         if (n == 0)
1157                 return fold;
1158
1159         fold = isl_qpolynomial_fold_cow(fold);
1160         if (!fold)
1161                 return NULL;
1162
1163         for (i = 0; i < fold->n; ++i) {
1164                 fold->qp[i] = isl_qpolynomial_substitute(fold->qp[i],
1165                                 type, first, n, subs);
1166                 if (!fold->qp[i])
1167                         goto error;
1168         }
1169
1170         return fold;
1171 error:
1172         isl_qpolynomial_fold_free(fold);
1173         return NULL;
1174 }