isl_range.c

   1 #include <isl_constraint.h>
   2 #include <isl_set.h>
   3 #include <isl_polynomial_private.h>
   4 #include <isl_morph.h>
   5
   6 struct range_data {
   7         int                     *signs;
   8         int                     sign;
   9         int                     test_monotonicity;
  10         int                     monotonicity;
  11         int                     exact;
  12         isl_qpolynomial         *qp;
  13         isl_qpolynomial         *poly;
  14         isl_pw_qpolynomial_fold *pwf;
  15         isl_pw_qpolynomial_fold *pwf_exact;
  16 };
  17
  18 static int propagate_on_domain(__isl_take isl_basic_set *bset,
  19         __isl_take isl_qpolynomial *poly, struct range_data *data);
  20
  21 /* Check whether the polynomial "poly" has sign "sign" over "bset",
  22  * i.e., if sign == 1, check that the lower bound on the polynomial
  23  * is non-negative and if sign == -1, check that the upper bound on
  24  * the polynomial is non-positive.
  25  */
  26 static int has_sign(__isl_keep isl_basic_set *bset,
  27         __isl_keep isl_qpolynomial *poly, int sign, int *signs)
  28 {
  29         struct range_data data_m;
  30         unsigned nvar;
  31         unsigned nparam;
  32         isl_dim *dim;
  33         isl_qpolynomial *opt;
  34         int r;
  35
  36         nparam = isl_basic_set_dim(bset, isl_dim_param);
  37         nvar = isl_basic_set_dim(bset, isl_dim_set);
  38
  39         bset = isl_basic_set_copy(bset);
  40         poly = isl_qpolynomial_copy(poly);
  41
  42         bset = isl_basic_set_move_dims(bset, isl_dim_set, 0,
  43                                         isl_dim_param, 0, nparam);
  44         poly = isl_qpolynomial_move_dims(poly, isl_dim_set, 0,
  45                                         isl_dim_param, 0, nparam);
  46
  47         dim = isl_qpolynomial_get_dim(poly);
  48         dim = isl_dim_drop(dim, isl_dim_set, 0, isl_dim_size(dim, isl_dim_set));
  49
  50         data_m.test_monotonicity = 0;
  51         data_m.signs = signs;
  52         data_m.pwf = isl_pw_qpolynomial_fold_zero(dim);
  53         data_m.sign = -sign;
  54         data_m.exact = 0;
  55         data_m.pwf_exact = NULL;
  56
  57         if (propagate_on_domain(bset, poly, &data_m) < 0)
  58                 goto error;
  59
  60         if (sign > 0)
  61                 opt = isl_pw_qpolynomial_fold_min(data_m.pwf);
  62         else
  63                 opt = isl_pw_qpolynomial_fold_max(data_m.pwf);
  64
  65         if (!opt)
  66                 r = -1;
  67         else if (isl_qpolynomial_is_nan(opt) ||
  68                  isl_qpolynomial_is_infty(opt) ||
  69                  isl_qpolynomial_is_neginfty(opt))
  70                 r = 0;
  71         else
  72                 r = sign * isl_qpolynomial_sgn(opt) >= 0;
  73
  74         isl_qpolynomial_free(opt);
  75
  76         return r;
  77 error:
  78         isl_pw_qpolynomial_fold_free(data_m.pwf);
  79         return -1;
  80 }
  81
  82 /* Return  1 if poly is monotonically increasing in the last set variable,
  83  *        -1 if poly is monotonically decreasing in the last set variable,
  84  *         0 if no conclusion,
  85  *        -2 on error.
  86  *
  87  * We simply check the sign of p(x+1)-p(x)
  88  */
  89 static int monotonicity(__isl_keep isl_basic_set *bset,
  90         __isl_keep isl_qpolynomial *poly, struct range_data *data)
  91 {
  92         isl_ctx *ctx;
  93         isl_dim *dim;
  94         isl_qpolynomial *sub = NULL;
  95         isl_qpolynomial *diff = NULL;
  96         int result = 0;
  97         int s;
  98         unsigned nvar;
  99
 100         ctx = isl_qpolynomial_get_ctx(poly);
 101         dim = isl_qpolynomial_get_dim(poly);
 102
 103         nvar = isl_basic_set_dim(bset, isl_dim_set);
 104
 105         sub = isl_qpolynomial_var(isl_dim_copy(dim), isl_dim_set, nvar - 1);
 106         sub = isl_qpolynomial_add(sub,
 107                 isl_qpolynomial_rat_cst(dim, ctx->one, ctx->one));
 108
 109         diff = isl_qpolynomial_substitute(isl_qpolynomial_copy(poly),
 110                         isl_dim_set, nvar - 1, 1, &sub);
 111         diff = isl_qpolynomial_sub(diff, isl_qpolynomial_copy(poly));
 112
 113         s = has_sign(bset, diff, 1, data->signs);
 114         if (s < 0)
 115                 goto error;
 116         if (s)
 117                 result = 1;
 118         else {
 119                 s = has_sign(bset, diff, -1, data->signs);
 120                 if (s < 0)
 121                         goto error;
 122                 if (s)
 123                         result = -1;
 124         }
 125
 126         isl_qpolynomial_free(diff);
 127         isl_qpolynomial_free(sub);
 128
 129         return result;
 130 error:
 131         isl_qpolynomial_free(diff);
 132         isl_qpolynomial_free(sub);
 133         return -2;
 134 }
 135
 136 static __isl_give isl_qpolynomial *bound2poly(__isl_take isl_constraint *bound,
 137         __isl_take isl_dim *dim, unsigned pos, int sign)
 138 {
 139         if (!bound) {
 140                 if (sign > 0)
 141                         return isl_qpolynomial_infty(dim);
 142                 else
 143                         return isl_qpolynomial_neginfty(dim);
 144         }
 145         isl_dim_free(dim);
 146         return isl_qpolynomial_from_constraint(bound, isl_dim_set, pos);
 147 }
 148
 149 static int bound_is_integer(__isl_take isl_constraint *bound, unsigned pos)
 150 {
 151         isl_int c;
 152         int is_int;
 153
 154         if (!bound)
 155                 return 1;
 156
 157         isl_int_init(c);
 158         isl_constraint_get_coefficient(bound, isl_dim_set, pos, &c);
 159         is_int = isl_int_is_one(c) || isl_int_is_negone(c);
 160         isl_int_clear(c);
 161
 162         return is_int;
 163 }
 164
 165 struct isl_fixed_sign_data {
 166         int             *signs;
 167         int             sign;
 168         isl_qpolynomial *poly;
 169 };
 170
 171 /* Add term "term" to data->poly if it has sign data->sign.
 172  * The sign is determined based on the signs of the parameters
 173  * and variables in data->signs.
 174  */
 175 static int collect_fixed_sign_terms(__isl_take isl_term *term, void *user)
 176 {
 177         struct isl_fixed_sign_data *data = (struct isl_fixed_sign_data *)user;
 178         isl_int n, d;
 179         int i;
 180         int sign;
 181         unsigned nparam;
 182         unsigned nvar;
 183
 184         if (!term)
 185                 return -1;
 186
 187         nparam = isl_term_dim(term, isl_dim_param);
 188         nvar = isl_term_dim(term, isl_dim_set);
 189
 190         isl_assert(isl_term_get_ctx(term), isl_term_dim(term, isl_dim_div) == 0,
 191                         return -1);
 192
 193         isl_int_init(n);
 194         isl_int_init(d);
 195
 196         isl_term_get_num(term, &n);
 197         isl_term_get_den(term, &d);
 198
 199         sign = isl_int_sgn(n);
 200         for (i = 0; i < nparam; ++i) {
 201                 if (data->signs[i] > 0)
 202                         continue;
 203                 if (isl_term_get_exp(term, isl_dim_param, i) % 2)
 204                         sign = -sign;
 205         }
 206         for (i = 0; i < nvar; ++i) {
 207                 if (data->signs[nparam + i] > 0)
 208                         continue;
 209                 if (isl_term_get_exp(term, isl_dim_set, i) % 2)
 210                         sign = -sign;
 211         }
 212
 213         if (sign == data->sign) {
 214                 isl_qpolynomial *t = isl_qpolynomial_from_term(term);
 215
 216                 data->poly = isl_qpolynomial_add(data->poly, t);
 217         } else
 218                 isl_term_free(term);
 219
 220         isl_int_clear(n);
 221         isl_int_clear(d);
 222
 223         return 0;
 224 }
 225
 226 /* Construct and return a polynomial that consists of the terms
 227  * in "poly" that have sign "sign".
 228  */
 229 static __isl_give isl_qpolynomial *fixed_sign_terms(
 230         __isl_keep isl_qpolynomial *poly, int *signs, int sign)
 231 {
 232         struct isl_fixed_sign_data data = { signs, sign };
 233         data.poly = isl_qpolynomial_zero(isl_qpolynomial_get_dim(poly));
 234
 235         if (isl_qpolynomial_foreach_term(poly, collect_fixed_sign_terms, &data) < 0)
 236                 goto error;
 237
 238         return data.poly;
 239 error:
 240         isl_qpolynomial_free(data.poly);
 241         return NULL;
 242 }
 243
 244 /* Helper function to add a guarder polynomial to either pwf_exact or pwf,
 245  * depending on whether the result has been determined to be exact.
 246  */
 247 static int add_guarded_poly(__isl_take isl_basic_set *bset,
 248         __isl_take isl_qpolynomial *poly, struct range_data *data)
 249 {
 250         enum isl_fold type = data->sign < 0 ? isl_fold_min : isl_fold_max;
 251         isl_set *set;
 252         isl_qpolynomial_fold *fold;
 253         isl_pw_qpolynomial_fold *pwf;
 254
 255         fold = isl_qpolynomial_fold_alloc(type, poly);
 256         set = isl_set_from_basic_set(bset);
 257         pwf = isl_pw_qpolynomial_fold_alloc(set, fold);
 258         if (data->exact)
 259                 data->pwf_exact = isl_pw_qpolynomial_fold_add(
 260                                                 data->pwf_exact, pwf);
 261         else
 262                 data->pwf = isl_pw_qpolynomial_fold_add(data->pwf, pwf);
 263
 264         return 0;
 265 }
 266
 267 /* Given a lower and upper bound on the final variable and constraints
 268  * on the remaining variables where these bounds are active,
 269  * eliminate the variable from data->poly based on these bounds.
 270  * If the polynomial has been determined to be monotonic
 271  * in the variable, then simply plug in the appropriate bound.
 272  * If the current polynomial is exact and if this bound is integer,
 273  * then the result is still exact.  In all other cases, the results
 274  * may be inexact.
 275  * Otherwise, plug in the largest bound (in absolute value) in
 276  * the positive terms (if an upper bound is wanted) or the negative terms
 277  * (if a lower bounded is wanted) and the other bound in the other terms.
 278  *
 279  * If all variables have been eliminated, then record the result.
 280  * Ohterwise, recurse on the next variable.
 281  */
 282 static int propagate_on_bound_pair(__isl_take isl_constraint *lower,
 283         __isl_take isl_constraint *upper, __isl_take isl_basic_set *bset,
 284         void *user)
 285 {
 286         struct range_data *data = (struct range_data *)user;
 287         int save_exact = data->exact;
 288         isl_qpolynomial *poly;
 289         int r;
 290         unsigned nvar;
 291
 292         nvar = isl_basic_set_dim(bset, isl_dim_set);
 293
 294         if (data->monotonicity) {
 295                 isl_qpolynomial *sub;
 296                 isl_dim *dim = isl_qpolynomial_get_dim(data->poly);
 297                 if (data->monotonicity * data->sign > 0) {
 298                         if (data->exact)
 299                                 data->exact = bound_is_integer(upper, nvar);
 300                         sub = bound2poly(upper, dim, nvar, 1);
 301                         isl_constraint_free(lower);
 302                 } else {
 303                         if (data->exact)
 304                                 data->exact = bound_is_integer(lower, nvar);
 305                         sub = bound2poly(lower, dim, nvar, -1);
 306                         isl_constraint_free(upper);
 307                 }
 308                 poly = isl_qpolynomial_copy(data->poly);
 309                 poly = isl_qpolynomial_substitute(poly, isl_dim_set, nvar, 1, &sub);
 310                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, nvar, 1);
 311
 312                 isl_qpolynomial_free(sub);
 313         } else {
 314                 isl_qpolynomial *l, *u;
 315                 isl_qpolynomial *pos, *neg;
 316                 isl_dim *dim = isl_qpolynomial_get_dim(data->poly);
 317                 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 318                 int sign = data->sign * data->signs[nparam + nvar];
 319
 320                 data->exact = 0;
 321
 322                 u = bound2poly(upper, isl_dim_copy(dim), nvar, 1);
 323                 l = bound2poly(lower, dim, nvar, -1);
 324
 325                 pos = fixed_sign_terms(data->poly, data->signs, sign);
 326                 neg = fixed_sign_terms(data->poly, data->signs, -sign);
 327
 328                 pos = isl_qpolynomial_substitute(pos, isl_dim_set, nvar, 1, &u);
 329                 neg = isl_qpolynomial_substitute(neg, isl_dim_set, nvar, 1, &l);
 330
 331                 poly = isl_qpolynomial_add(pos, neg);
 332                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, nvar, 1);
 333
 334                 isl_qpolynomial_free(u);
 335                 isl_qpolynomial_free(l);
 336         }
 337
 338         if (isl_basic_set_dim(bset, isl_dim_set) == 0)
 339                 r = add_guarded_poly(bset, poly, data);
 340         else
 341                 r = propagate_on_domain(bset, poly, data);
 342
 343         data->exact = save_exact;
 344
 345         return r;
 346 }
 347
 348 /* Recursively perform range propagation on the polynomial "poly"
 349  * defined over the basic set "bset" and collect the results in "data".
 350  */
 351 static int propagate_on_domain(__isl_take isl_basic_set *bset,
 352         __isl_take isl_qpolynomial *poly, struct range_data *data)
 353 {
 354         isl_qpolynomial *save_poly = data->poly;
 355         int save_monotonicity = data->monotonicity;
 356         unsigned d;
 357
 358         if (!bset || !poly)
 359                 goto error;
 360
 361         d = isl_basic_set_dim(bset, isl_dim_set);
 362         isl_assert(bset->ctx, d >= 1, goto error);
 363
 364         if (isl_qpolynomial_is_cst(poly, NULL, NULL)) {
 365                 bset = isl_basic_set_project_out(bset, isl_dim_set, 0, d);
 366                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, 0, d);
 367                 return add_guarded_poly(bset, poly, data);
 368         }
 369
 370         if (data->test_monotonicity)
 371                 data->monotonicity = monotonicity(bset, poly, data);
 372         else
 373                 data->monotonicity = 0;
 374         if (data->monotonicity < -1)
 375                 goto error;
 376
 377         data->poly = poly;
 378         if (isl_basic_set_foreach_bound_pair(bset, isl_dim_set, d - 1,
 379                                             &propagate_on_bound_pair, data) < 0)
 380                 goto error;
 381
 382         isl_basic_set_free(bset);
 383         isl_qpolynomial_free(poly);
 384         data->monotonicity = save_monotonicity;
 385         data->poly = save_poly;
 386
 387         return 0;
 388 error:
 389         isl_basic_set_free(bset);
 390         isl_qpolynomial_free(poly);
 391         data->monotonicity = save_monotonicity;
 392         data->poly = save_poly;
 393         return -1;
 394 }
 395
 396 static int basic_guarded_poly_bound(__isl_take isl_basic_set *bset, void *user)
 397 {
 398         struct range_data *data = (struct range_data *)user;
 399         unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 400         unsigned dim = isl_basic_set_dim(bset, isl_dim_set);
 401         int r;
 402
 403         data->signs = NULL;
 404
 405         data->signs = isl_alloc_array(bset->ctx, int,
 406                                         isl_basic_set_dim(bset, isl_dim_all));
 407
 408         if (isl_basic_set_dims_get_sign(bset, isl_dim_set, 0, dim,
 409                                         data->signs + nparam) < 0)
 410                 goto error;
 411         if (isl_basic_set_dims_get_sign(bset, isl_dim_param, 0, nparam,
 412                                         data->signs) < 0)
 413                 goto error;
 414
 415         r = propagate_on_domain(bset, isl_qpolynomial_copy(data->poly), data);
 416
 417         free(data->signs);
 418
 419         return r;
 420 error:
 421         free(data->signs);
 422         isl_basic_set_free(bset);
 423         return -1;
 424 }
 425
 426 static int compressed_guarded_poly_bound(__isl_take isl_basic_set *bset,
 427         __isl_take isl_qpolynomial *poly, void *user)
 428 {
 429         struct range_data *data = (struct range_data *)user;
 430         unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 431         unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
 432         isl_set *set;
 433
 434         if (!bset)
 435                 goto error;
 436
 437         if (nvar == 0)
 438                 return add_guarded_poly(bset, poly, data);
 439
 440         set = isl_set_from_basic_set(bset);
 441         set = isl_set_split_dims(set, isl_dim_param, 0, nparam);
 442         set = isl_set_split_dims(set, isl_dim_set, 0, nvar);
 443
 444         data->poly = poly;
 445
 446         if (isl_set_foreach_basic_set(set, &basic_guarded_poly_bound, data) < 0)
 447                 goto error;
 448
 449         isl_set_free(set);
 450         isl_qpolynomial_free(poly);
 451
 452         return 0;
 453 error:
 454         isl_set_free(set);
 455         isl_qpolynomial_free(poly);
 456         return -1;
 457 }
 458
 459 static int guarded_poly_bound(__isl_take isl_basic_set *bset,
 460         __isl_take isl_qpolynomial *poly, void *user)
 461 {
 462         struct range_data *data = (struct range_data *)user;
 463         isl_pw_qpolynomial_fold *top_pwf;
 464         isl_pw_qpolynomial_fold *top_pwf_exact;
 465         isl_dim *dim;
 466         isl_morph *morph, *morph2;
 467         unsigned orig_nvar;
 468         int r;
 469
 470         bset = isl_basic_set_detect_equalities(bset);
 471
 472         if (!bset)
 473                 goto error;
 474
 475         if (bset->n_eq == 0)
 476                 return compressed_guarded_poly_bound(bset, poly, user);
 477
 478         orig_nvar = isl_basic_set_dim(bset, isl_dim_set);
 479
 480         morph = isl_basic_set_variable_compression(bset, isl_dim_param);
 481         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
 482
 483         morph2 = isl_basic_set_parameter_compression(bset);
 484         bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
 485
 486         morph = isl_morph_compose(morph2, morph);
 487
 488         morph2 = isl_basic_set_variable_compression(bset, isl_dim_set);
 489         bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
 490
 491         morph2 = isl_morph_compose(morph2, isl_morph_copy(morph));
 492         poly = isl_qpolynomial_morph(poly, morph2);
 493
 494         dim = isl_morph_get_ran_dim(morph);
 495         dim = isl_dim_drop(dim, isl_dim_set, 0, isl_dim_size(dim, isl_dim_set));
 496
 497         top_pwf = data->pwf;
 498         top_pwf_exact = data->pwf_exact;
 499
 500         data->pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
 501         data->pwf_exact = isl_pw_qpolynomial_fold_zero(dim);
 502
 503         r = compressed_guarded_poly_bound(bset, poly, user);
 504
 505         morph = isl_morph_remove_dom_dims(morph, isl_dim_set, 0, orig_nvar);
 506         morph = isl_morph_remove_ran_dims(morph, isl_dim_set, 0, orig_nvar);
 507         morph = isl_morph_inverse(morph);
 508
 509         data->pwf = isl_pw_qpolynomial_fold_morph(data->pwf,
 510                                                         isl_morph_copy(morph));
 511         data->pwf_exact = isl_pw_qpolynomial_fold_morph(data->pwf_exact, morph);
 512
 513         data->pwf = isl_pw_qpolynomial_fold_add(top_pwf, data->pwf);
 514         data->pwf_exact = isl_pw_qpolynomial_fold_add(top_pwf_exact,
 515                                                         data->pwf_exact);
 516
 517         return r;
 518 error:
 519         isl_basic_set_free(bset);
 520         isl_qpolynomial_free(poly);
 521         return -1;
 522 }
 523
 524 static int basic_guarded_bound(__isl_take isl_basic_set *bset, void *user)
 525 {
 526         struct range_data *data = (struct range_data *)user;
 527         int r;
 528
 529         r = isl_qpolynomial_as_polynomial_on_domain(data->qp, bset,
 530                                                     &guarded_poly_bound, user);
 531         isl_basic_set_free(bset);
 532         return r;
 533 }
 534
 535 static int guarded_bound(__isl_take isl_set *set,
 536         __isl_take isl_qpolynomial *qp, void *user)
 537 {
 538         struct range_data *data = (struct range_data *)user;
 539
 540         if (!set || !qp)
 541                 goto error;
 542
 543         set = isl_set_make_disjoint(set);
 544
 545         data->qp = qp;
 546
 547         if (isl_set_foreach_basic_set(set, &basic_guarded_bound, data) < 0)
 548                 goto error;
 549
 550         isl_set_free(set);
 551         isl_qpolynomial_free(qp);
 552
 553         return 0;
 554 error:
 555         isl_set_free(set);
 556         isl_qpolynomial_free(qp);
 557         return -1;
 558 }
 559
 560 /* Compute a lower or upper bound (depending on "type") on the given
 561  * piecewise step-polynomial using range propagation.
 562  */
 563 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound_range(
 564         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type, int *exact)
 565 {
 566         isl_dim *dim;
 567         isl_pw_qpolynomial_fold *pwf;
 568         unsigned nvar;
 569         unsigned nparam;
 570         struct range_data data;
 571         int covers;
 572
 573         if (!pwqp)
 574                 return NULL;
 575
 576         dim = isl_pw_qpolynomial_get_dim(pwqp);
 577         nvar = isl_dim_size(dim, isl_dim_set);
 578
 579         if (isl_pw_qpolynomial_is_zero(pwqp)) {
 580                 isl_pw_qpolynomial_free(pwqp);
 581                 dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 582                 return isl_pw_qpolynomial_fold_zero(dim);
 583         }
 584
 585         if (nvar == 0) {
 586                 isl_dim_free(dim);
 587                 return isl_pw_qpolynomial_fold_from_pw_qpolynomial(type, pwqp);
 588         }
 589
 590         dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 591
 592         nparam = isl_dim_size(dim, isl_dim_param);
 593         data.pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
 594         data.pwf_exact = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
 595         if (type == isl_fold_min)
 596                 data.sign = -1;
 597         else
 598                 data.sign = 1;
 599         data.test_monotonicity = 1;
 600         data.exact = !!exact;
 601
 602         if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp, guarded_bound, &data))
 603                 goto error;
 604
 605         covers = isl_pw_qpolynomial_fold_covers(data.pwf_exact, data.pwf);
 606         if (covers < 0)
 607                 goto error;
 608
 609         if (exact)
 610                 *exact = covers;
 611
 612         isl_dim_free(dim);
 613         isl_pw_qpolynomial_free(pwqp);
 614
 615         if (covers) {
 616                 isl_pw_qpolynomial_fold_free(data.pwf);
 617                 return data.pwf_exact;
 618         }
 619
 620         data.pwf = isl_pw_qpolynomial_fold_add(data.pwf, data.pwf_exact);
 621
 622         return data.pwf;
 623 error:
 624         isl_pw_qpolynomial_fold_free(data.pwf);
 625         isl_dim_free(dim);
 626         isl_pw_qpolynomial_free(pwqp);
 627         return NULL;
 628 }