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_polynomial_private.h>
  12 #include <isl_point_private.h>
  13 #include <isl_dim.h>
  14 #include <isl_map_private.h>
  15 #include <isl_lp.h>
  16 #include <isl_seq.h>
  17
  18 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
  19         enum isl_fold type, __isl_take isl_dim *dim, int n)
  20 {
  21         isl_qpolynomial_fold *fold;
  22
  23         if (!dim)
  24                 goto error;
  25
  26         isl_assert(dim->ctx, n >= 0, goto error);
  27         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
  28                         sizeof(struct isl_qpolynomial_fold) +
  29                         (n - 1) * sizeof(struct isl_qpolynomial *));
  30         if (!fold)
  31                 goto error;
  32
  33         fold->ref = 1;
  34         fold->size = n;
  35         fold->n = 0;
  36         fold->type = type;
  37         fold->dim = dim;
  38
  39         return fold;
  40 error:
  41         isl_dim_free(dim);
  42         return NULL;
  43 }
  44
  45 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
  46         enum isl_dim_type type, unsigned first, unsigned n)
  47 {
  48         int i;
  49
  50         if (!fold)
  51                 return -1;
  52         if (fold->n == 0 || n == 0)
  53                 return 0;
  54
  55         for (i = 0; i < fold->n; ++i) {
  56                 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
  57                                                             type, first, n);
  58                 if (involves < 0 || involves)
  59                         return involves;
  60         }
  61         return 0;
  62 }
  63
  64 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
  65         __isl_take isl_qpolynomial_fold *fold,
  66         enum isl_dim_type type, unsigned first, unsigned n)
  67 {
  68         int i;
  69
  70         if (!fold)
  71                 return NULL;
  72         if (n == 0)
  73                 return fold;
  74
  75         fold = isl_qpolynomial_fold_cow(fold);
  76         if (!fold)
  77                 return NULL;
  78         fold->dim = isl_dim_drop(fold->dim, type, first, n);
  79         if (!fold->dim)
  80                 goto error;
  81
  82         for (i = 0; i < fold->n; ++i) {
  83                 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
  84                                                             type, first, n);
  85                 if (!fold->qp[i])
  86                         goto error;
  87         }
  88
  89         return fold;
  90 error:
  91         isl_qpolynomial_fold_free(fold);
  92         return NULL;
  93 }
  94
  95 static int isl_qpolynomial_cst_sign(__isl_keep isl_qpolynomial *qp)
  96 {
  97         struct isl_upoly_cst *cst;
  98
  99         cst = isl_upoly_as_cst(qp->upoly);
 100         if (!cst)
 101                 return 0;
 102
 103         return isl_int_sgn(cst->n) < 0 ? -1 : 1;
 104 }
 105
 106 static int isl_qpolynomial_aff_sign(__isl_keep isl_set *set,
 107         __isl_keep isl_qpolynomial *qp)
 108 {
 109         enum isl_lp_result res;
 110         isl_vec *aff;
 111         isl_int opt;
 112         int sgn = 0;
 113
 114         aff = isl_qpolynomial_extract_affine(qp);
 115         if (!aff)
 116                 return 0;
 117
 118         isl_int_init(opt);
 119
 120         res = isl_set_solve_lp(set, 0, aff->el + 1, aff->el[0],
 121                                 &opt, NULL, NULL);
 122         if (res == isl_lp_error)
 123                 goto done;
 124         if (res == isl_lp_empty ||
 125             (res == isl_lp_ok && !isl_int_is_neg(opt))) {
 126                 sgn = 1;
 127                 goto done;
 128         }
 129
 130         res = isl_set_solve_lp(set, 1, aff->el + 1, aff->el[0],
 131                                 &opt, NULL, NULL);
 132         if (res == isl_lp_ok && !isl_int_is_pos(opt))
 133                 sgn = -1;
 134
 135 done:
 136         isl_int_clear(opt);
 137         isl_vec_free(aff);
 138         return sgn;
 139 }
 140
 141 /* Determine, if possible, the sign of the quasipolynomial "qp" on
 142  * the domain "set".
 143  *
 144  * If qp is a constant, then the problem is trivial.
 145  * If qp is linear, then we check if the minimum of the corresponding
 146  * affine constraint is non-negative or if the maximum is non-positive.
 147  *
 148  * Otherwise, we check if the outermost variable "v" has a lower bound "l"
 149  * in "set".  If so, we write qp(v,v') as
 150  *
 151  *      q(v,v') * (v - l) + r(v')
 152  *
 153  * if q(v,v') and r(v') have the same known sign, then the original
 154  * quasipolynomial has the same sign as well.
 155  *
 156  * Return
 157  *      -1 if qp <= 0
 158  *       1 if qp >= 0
 159  *       0 if unknown
 160  */
 161 static int isl_qpolynomial_sign(__isl_keep isl_set *set,
 162         __isl_keep isl_qpolynomial *qp)
 163 {
 164         int d;
 165         int i;
 166         int is;
 167         struct isl_upoly_rec *rec;
 168         isl_vec *v;
 169         isl_int l;
 170         enum isl_lp_result res;
 171         int sgn = 0;
 172
 173         is = isl_qpolynomial_is_cst(qp, NULL, NULL);
 174         if (is < 0)
 175                 return 0;
 176         if (is)
 177                 return isl_qpolynomial_cst_sign(qp);
 178
 179         is = isl_qpolynomial_is_affine(qp);
 180         if (is < 0)
 181                 return 0;
 182         if (is)
 183                 return isl_qpolynomial_aff_sign(set, qp);
 184
 185         if (qp->div->n_row > 0)
 186                 return 0;
 187
 188         rec = isl_upoly_as_rec(qp->upoly);
 189         if (!rec)
 190                 return 0;
 191
 192         d = isl_dim_total(qp->dim);
 193         v = isl_vec_alloc(set->ctx, 2 + d);
 194         if (!v)
 195                 return 0;
 196
 197         isl_seq_clr(v->el + 1, 1 + d);
 198         isl_int_set_si(v->el[0], 1);
 199         isl_int_set_si(v->el[2 + qp->upoly->var], 1);
 200
 201         isl_int_init(l);
 202
 203         res = isl_set_solve_lp(set, 0, v->el + 1, v->el[0], &l, NULL, NULL);
 204         if (res == isl_lp_ok) {
 205                 isl_qpolynomial *min;
 206                 isl_qpolynomial *base;
 207                 isl_qpolynomial *r, *q;
 208                 isl_qpolynomial *t;
 209
 210                 min = isl_qpolynomial_cst(isl_dim_copy(qp->dim), l);
 211                 base = isl_qpolynomial_pow(isl_dim_copy(qp->dim),
 212                                                 qp->upoly->var, 1);
 213
 214                 r = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
 215                                           isl_upoly_copy(rec->p[rec->n - 1]));
 216                 q = isl_qpolynomial_copy(r);
 217
 218                 for (i = rec->n - 2; i >= 0; --i) {
 219                         r = isl_qpolynomial_mul(r, isl_qpolynomial_copy(min));
 220                         t = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
 221                                                   isl_upoly_copy(rec->p[i]));
 222                         r = isl_qpolynomial_add(r, t);
 223                         if (i == 0)
 224                                 break;
 225                         q = isl_qpolynomial_mul(q, isl_qpolynomial_copy(base));
 226                         q = isl_qpolynomial_add(q, isl_qpolynomial_copy(r));
 227                 }
 228
 229                 if (isl_qpolynomial_is_zero(q))
 230                         sgn = isl_qpolynomial_sign(set, r);
 231                 else if (isl_qpolynomial_is_zero(r))
 232                         sgn = isl_qpolynomial_sign(set, q);
 233                 else {
 234                         int sgn_q, sgn_r;
 235                         sgn_r = isl_qpolynomial_sign(set, r);
 236                         sgn_q = isl_qpolynomial_sign(set, q);
 237                         if (sgn_r == sgn_q)
 238                                 sgn = sgn_r;
 239                 }
 240
 241                 isl_qpolynomial_free(min);
 242                 isl_qpolynomial_free(base);
 243                 isl_qpolynomial_free(q);
 244                 isl_qpolynomial_free(r);
 245         }
 246
 247         isl_int_clear(l);
 248
 249         isl_vec_free(v);
 250
 251         return sgn;
 252 }
 253
 254 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold_on_domain(
 255         __isl_keep isl_set *set,
 256         __isl_take isl_qpolynomial_fold *fold1,
 257         __isl_take isl_qpolynomial_fold *fold2)
 258 {
 259         int i, j;
 260         int n1;
 261         struct isl_qpolynomial_fold *res = NULL;
 262         int better;
 263
 264         if (!fold1 || !fold2)
 265                 goto error;
 266
 267         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
 268         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
 269                         goto error);
 270
 271         better = fold1->type == isl_fold_max ? -1 : 1;
 272
 273         if (isl_qpolynomial_fold_is_empty(fold1)) {
 274                 isl_qpolynomial_fold_free(fold1);
 275                 return fold2;
 276         }
 277
 278         if (isl_qpolynomial_fold_is_empty(fold2)) {
 279                 isl_qpolynomial_fold_free(fold2);
 280                 return fold1;
 281         }
 282
 283         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
 284                                         fold1->n + fold2->n);
 285         if (!res)
 286                 goto error;
 287
 288         for (i = 0; i < fold1->n; ++i) {
 289                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
 290                 if (!res->qp[res->n])
 291                         goto error;
 292                 res->n++;
 293         }
 294         n1 = res->n;
 295
 296         for (i = 0; i < fold2->n; ++i) {
 297                 for (j = n1 - 1; j >= 0; --j) {
 298                         isl_qpolynomial *d;
 299                         int sgn;
 300                         d = isl_qpolynomial_sub(
 301                                 isl_qpolynomial_copy(res->qp[j]),
 302                                 isl_qpolynomial_copy(fold2->qp[i]));
 303                         sgn = isl_qpolynomial_sign(set, d);
 304                         isl_qpolynomial_free(d);
 305                         if (sgn == 0)
 306                                 continue;
 307                         if (sgn != better)
 308                                 break;
 309                         isl_qpolynomial_free(res->qp[j]);
 310                         if (j != n1 - 1)
 311                                 res->qp[j] = res->qp[n1 - 1];
 312                         n1--;
 313                         if (n1 != res->n - 1)
 314                                 res->qp[n1] = res->qp[res->n - 1];
 315                         res->n--;
 316                 }
 317                 if (j >= 0)
 318                         continue;
 319                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
 320                 if (!res->qp[res->n])
 321                         goto error;
 322                 res->n++;
 323         }
 324
 325         isl_qpolynomial_fold_free(fold1);
 326         isl_qpolynomial_fold_free(fold2);
 327
 328         return res;
 329 error:
 330         isl_qpolynomial_fold_free(res);
 331         isl_qpolynomial_fold_free(fold1);
 332         isl_qpolynomial_fold_free(fold2);
 333         return NULL;
 334 }
 335
 336 #undef PW
 337 #define PW isl_pw_qpolynomial_fold
 338 #undef EL
 339 #define EL isl_qpolynomial_fold
 340 #undef IS_ZERO
 341 #define IS_ZERO is_empty
 342 #undef FIELD
 343 #define FIELD fold
 344 #undef ADD
 345 #define ADD(d,e1,e2)    isl_qpolynomial_fold_fold_on_domain(d,e1,e2)
 346
 347 #include <isl_pw_templ.c>
 348
 349 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
 350         __isl_take isl_dim *dim)
 351 {
 352         return qpolynomial_fold_alloc(type, dim, 0);
 353 }
 354
 355 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
 356         enum isl_fold type, __isl_take isl_qpolynomial *qp)
 357 {
 358         isl_qpolynomial_fold *fold;
 359
 360         if (!qp)
 361                 return NULL;
 362
 363         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
 364         if (!fold)
 365                 goto error;
 366
 367         fold->qp[0] = qp;
 368         fold->n++;
 369
 370         return fold;
 371 error:
 372         isl_qpolynomial_fold_free(fold);
 373         isl_qpolynomial_free(qp);
 374         return NULL;
 375 }
 376
 377 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
 378         __isl_keep isl_qpolynomial_fold *fold)
 379 {
 380         if (!fold)
 381                 return NULL;
 382
 383         fold->ref++;
 384         return fold;
 385 }
 386
 387 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
 388         __isl_keep isl_qpolynomial_fold *fold)
 389 {
 390         int i;
 391         isl_qpolynomial_fold *dup;
 392
 393         if (!fold)
 394                 return NULL;
 395         dup = qpolynomial_fold_alloc(fold->type,
 396                                         isl_dim_copy(fold->dim), fold->n);
 397         if (!dup)
 398                 return NULL;
 399
 400         for (i = 0; i < fold->n; ++i) {
 401                 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
 402                 if (!dup->qp[i])
 403                         goto error;
 404         }
 405
 406         return dup;
 407 error:
 408         isl_qpolynomial_fold_free(dup);
 409         return NULL;
 410 }
 411
 412 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
 413         __isl_take isl_qpolynomial_fold *fold)
 414 {
 415         if (!fold)
 416                 return NULL;
 417
 418         if (fold->ref == 1)
 419                 return fold;
 420         fold->ref--;
 421         return isl_qpolynomial_fold_dup(fold);
 422 }
 423
 424 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
 425 {
 426         int i;
 427
 428         if (!fold)
 429                 return;
 430         if (--fold->ref > 0)
 431                 return;
 432
 433         for (i = 0; i < fold->n; ++i)
 434                 isl_qpolynomial_free(fold->qp[i]);
 435         isl_dim_free(fold->dim);
 436         free(fold);
 437 }
 438
 439 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
 440 {
 441         if (!fold)
 442                 return -1;
 443
 444         return fold->n == 0;
 445 }
 446
 447 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
 448         __isl_take isl_qpolynomial_fold *fold1,
 449         __isl_take isl_qpolynomial_fold *fold2)
 450 {
 451         int i;
 452         struct isl_qpolynomial_fold *res = NULL;
 453
 454         if (!fold1 || !fold2)
 455                 goto error;
 456
 457         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
 458         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
 459                         goto error);
 460
 461         if (isl_qpolynomial_fold_is_empty(fold1)) {
 462                 isl_qpolynomial_fold_free(fold1);
 463                 return fold2;
 464         }
 465
 466         if (isl_qpolynomial_fold_is_empty(fold2)) {
 467                 isl_qpolynomial_fold_free(fold2);
 468                 return fold1;
 469         }
 470
 471         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
 472                                         fold1->n + fold2->n);
 473         if (!res)
 474                 goto error;
 475
 476         for (i = 0; i < fold1->n; ++i) {
 477                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
 478                 if (!res->qp[res->n])
 479                         goto error;
 480                 res->n++;
 481         }
 482
 483         for (i = 0; i < fold2->n; ++i) {
 484                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
 485                 if (!res->qp[res->n])
 486                         goto error;
 487                 res->n++;
 488         }
 489
 490         isl_qpolynomial_fold_free(fold1);
 491         isl_qpolynomial_fold_free(fold2);
 492
 493         return res;
 494 error:
 495         isl_qpolynomial_fold_free(res);
 496         isl_qpolynomial_fold_free(fold1);
 497         isl_qpolynomial_fold_free(fold2);
 498         return NULL;
 499 }
 500
 501 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
 502         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
 503 {
 504         int i;
 505         isl_pw_qpolynomial_fold *pwf;
 506
 507         if (!pwqp)
 508                 return NULL;
 509
 510         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
 511
 512         for (i = 0; i < pwqp->n; ++i)
 513                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
 514                         isl_set_copy(pwqp->p[i].set),
 515                         isl_qpolynomial_fold_alloc(type,
 516                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
 517
 518         isl_pw_qpolynomial_free(pwqp);
 519
 520         return pwf;
 521 }
 522
 523 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
 524         __isl_keep isl_qpolynomial_fold *fold2)
 525 {
 526         int i;
 527
 528         if (!fold1 || !fold2)
 529                 return -1;
 530
 531         if (fold1->n != fold2->n)
 532                 return 0;
 533
 534         /* We probably want to sort the qps first... */
 535         for (i = 0; i < fold1->n; ++i) {
 536                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
 537                 if (eq < 0 || !eq)
 538                         return eq;
 539         }
 540
 541         return 1;
 542 }
 543
 544 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
 545         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
 546 {
 547         isl_qpolynomial *qp;
 548
 549         if (!fold || !pnt)
 550                 goto error;
 551         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
 552         isl_assert(pnt->dim->ctx,
 553                 fold->type == isl_fold_max || fold->type == isl_fold_min,
 554                 goto error);
 555
 556         if (fold->n == 0)
 557                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
 558         else {
 559                 int i;
 560                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
 561                                                 isl_point_copy(pnt));
 562                 for (i = 1; i < fold->n; ++i) {
 563                         isl_qpolynomial *qp_i;
 564                         qp_i = isl_qpolynomial_eval(
 565                                             isl_qpolynomial_copy(fold->qp[i]),
 566                                             isl_point_copy(pnt));
 567                         if (fold->type == isl_fold_max)
 568                                 qp = isl_qpolynomial_max_cst(qp, qp_i);
 569                         else
 570                                 qp = isl_qpolynomial_min_cst(qp, qp_i);
 571                 }
 572         }
 573         isl_qpolynomial_fold_free(fold);
 574         isl_point_free(pnt);
 575
 576         return qp;
 577 error:
 578         isl_qpolynomial_fold_free(fold);
 579         isl_point_free(pnt);
 580         return NULL;
 581 }
 582
 583 size_t isl_pw_qpolynomial_fold_size(__isl_keep isl_pw_qpolynomial_fold *pwf)
 584 {
 585         int i;
 586         size_t n = 0;
 587
 588         for (i = 0; i < pwf->n; ++i)
 589                 n += pwf->p[i].fold->n;
 590
 591         return n;
 592 }
 593
 594 __isl_give isl_qpolynomial *isl_qpolynomial_fold_opt_on_domain(
 595         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *set, int max)
 596 {
 597         int i;
 598         isl_qpolynomial *opt;
 599
 600         if (!set || !fold)
 601                 goto error;
 602
 603         if (fold->n == 0) {
 604                 isl_dim *dim = isl_dim_copy(fold->dim);
 605                 isl_set_free(set);
 606                 isl_qpolynomial_fold_free(fold);
 607                 return isl_qpolynomial_zero(dim);
 608         }
 609
 610         opt = isl_qpolynomial_opt_on_domain(isl_qpolynomial_copy(fold->qp[i]),
 611                                                 isl_set_copy(set), max);
 612         for (i = 1; i < fold->n; ++i) {
 613                 isl_qpolynomial *opt_i;
 614                 opt_i = isl_qpolynomial_opt_on_domain(
 615                                 isl_qpolynomial_copy(fold->qp[i]),
 616                                 isl_set_copy(set), max);
 617                 if (max)
 618                         opt = isl_qpolynomial_max_cst(opt, opt_i);
 619                 else
 620                         opt = isl_qpolynomial_min_cst(opt, opt_i);
 621         }
 622
 623         isl_set_free(set);
 624         isl_qpolynomial_fold_free(fold);
 625
 626         return opt;
 627 error:
 628         isl_set_free(set);
 629         isl_qpolynomial_fold_free(fold);
 630         return NULL;
 631 }