isl_range.c

   1 #include <isl_ctx_private.h>
   2 #include <isl/constraint.h>
   3 #include <isl/set.h>
   4 #include <isl_polynomial_private.h>
   5 #include <isl_morph.h>
   6 #include <isl_range.h>
   7
   8 struct range_data {
   9         struct isl_bound        *bound;
  10         int                     *signs;
  11         int                     sign;
  12         int                     test_monotonicity;
  13         int                     monotonicity;
  14         int                     tight;
  15         isl_qpolynomial         *poly;
  16         isl_pw_qpolynomial_fold *pwf;
  17         isl_pw_qpolynomial_fold *pwf_tight;
  18 };
  19
  20 static int propagate_on_domain(__isl_take isl_basic_set *bset,
  21         __isl_take isl_qpolynomial *poly, struct range_data *data);
  22
  23 /* Check whether the polynomial "poly" has sign "sign" over "bset",
  24  * i.e., if sign == 1, check that the lower bound on the polynomial
  25  * is non-negative and if sign == -1, check that the upper bound on
  26  * the polynomial is non-positive.
  27  */
  28 static int has_sign(__isl_keep isl_basic_set *bset,
  29         __isl_keep isl_qpolynomial *poly, int sign, int *signs)
  30 {
  31         struct range_data data_m;
  32         unsigned nvar;
  33         unsigned nparam;
  34         isl_space *dim;
  35         isl_qpolynomial *opt;
  36         int r;
  37         enum isl_fold type;
  38
  39         nparam = isl_basic_set_dim(bset, isl_dim_param);
  40         nvar = isl_basic_set_dim(bset, isl_dim_set);
  41
  42         bset = isl_basic_set_copy(bset);
  43         poly = isl_qpolynomial_copy(poly);
  44
  45         bset = isl_basic_set_move_dims(bset, isl_dim_set, 0,
  46                                         isl_dim_param, 0, nparam);
  47         poly = isl_qpolynomial_move_dims(poly, isl_dim_in, 0,
  48                                         isl_dim_param, 0, nparam);
  49
  50         dim = isl_qpolynomial_get_space(poly);
  51         dim = isl_space_params(dim);
  52         dim = isl_space_from_domain(dim);
  53         dim = isl_space_add_dims(dim, isl_dim_out, 1);
  54
  55         data_m.test_monotonicity = 0;
  56         data_m.signs = signs;
  57         data_m.sign = -sign;
  58         type = data_m.sign < 0 ? isl_fold_min : isl_fold_max;
  59         data_m.pwf = isl_pw_qpolynomial_fold_zero(dim, type);
  60         data_m.tight = 0;
  61         data_m.pwf_tight = NULL;
  62
  63         if (propagate_on_domain(bset, poly, &data_m) < 0)
  64                 goto error;
  65
  66         if (sign > 0)
  67                 opt = isl_pw_qpolynomial_fold_min(data_m.pwf);
  68         else
  69                 opt = isl_pw_qpolynomial_fold_max(data_m.pwf);
  70
  71         if (!opt)
  72                 r = -1;
  73         else if (isl_qpolynomial_is_nan(opt) ||
  74                  isl_qpolynomial_is_infty(opt) ||
  75                  isl_qpolynomial_is_neginfty(opt))
  76                 r = 0;
  77         else
  78                 r = sign * isl_qpolynomial_sgn(opt) >= 0;
  79
  80         isl_qpolynomial_free(opt);
  81
  82         return r;
  83 error:
  84         isl_pw_qpolynomial_fold_free(data_m.pwf);
  85         return -1;
  86 }
  87
  88 /* Return  1 if poly is monotonically increasing in the last set variable,
  89  *        -1 if poly is monotonically decreasing in the last set variable,
  90  *         0 if no conclusion,
  91  *        -2 on error.
  92  *
  93  * We simply check the sign of p(x+1)-p(x)
  94  */
  95 static int monotonicity(__isl_keep isl_basic_set *bset,
  96         __isl_keep isl_qpolynomial *poly, struct range_data *data)
  97 {
  98         isl_ctx *ctx;
  99         isl_space *dim;
 100         isl_qpolynomial *sub = NULL;
 101         isl_qpolynomial *diff = NULL;
 102         int result = 0;
 103         int s;
 104         unsigned nvar;
 105
 106         ctx = isl_qpolynomial_get_ctx(poly);
 107         dim = isl_qpolynomial_get_domain_space(poly);
 108
 109         nvar = isl_basic_set_dim(bset, isl_dim_set);
 110
 111         sub = isl_qpolynomial_var_on_domain(isl_space_copy(dim), isl_dim_set, nvar - 1);
 112         sub = isl_qpolynomial_add(sub,
 113                 isl_qpolynomial_rat_cst_on_domain(dim, ctx->one, ctx->one));
 114
 115         diff = isl_qpolynomial_substitute(isl_qpolynomial_copy(poly),
 116                         isl_dim_in, nvar - 1, 1, &sub);
 117         diff = isl_qpolynomial_sub(diff, isl_qpolynomial_copy(poly));
 118
 119         s = has_sign(bset, diff, 1, data->signs);
 120         if (s < 0)
 121                 goto error;
 122         if (s)
 123                 result = 1;
 124         else {
 125                 s = has_sign(bset, diff, -1, data->signs);
 126                 if (s < 0)
 127                         goto error;
 128                 if (s)
 129                         result = -1;
 130         }
 131
 132         isl_qpolynomial_free(diff);
 133         isl_qpolynomial_free(sub);
 134
 135         return result;
 136 error:
 137         isl_qpolynomial_free(diff);
 138         isl_qpolynomial_free(sub);
 139         return -2;
 140 }
 141
 142 static __isl_give isl_qpolynomial *bound2poly(__isl_take isl_constraint *bound,
 143         __isl_take isl_space *dim, unsigned pos, int sign)
 144 {
 145         if (!bound) {
 146                 if (sign > 0)
 147                         return isl_qpolynomial_infty_on_domain(dim);
 148                 else
 149                         return isl_qpolynomial_neginfty_on_domain(dim);
 150         }
 151         isl_space_free(dim);
 152         return isl_qpolynomial_from_constraint(bound, isl_dim_set, pos);
 153 }
 154
 155 static int bound_is_integer(__isl_take isl_constraint *bound, unsigned pos)
 156 {
 157         isl_int c;
 158         int is_int;
 159
 160         if (!bound)
 161                 return 1;
 162
 163         isl_int_init(c);
 164         isl_constraint_get_coefficient(bound, isl_dim_set, pos, &c);
 165         is_int = isl_int_is_one(c) || isl_int_is_negone(c);
 166         isl_int_clear(c);
 167
 168         return is_int;
 169 }
 170
 171 struct isl_fixed_sign_data {
 172         int             *signs;
 173         int             sign;
 174         isl_qpolynomial *poly;
 175 };
 176
 177 /* Add term "term" to data->poly if it has sign data->sign.
 178  * The sign is determined based on the signs of the parameters
 179  * and variables in data->signs.  The integer divisions, if
 180  * any, are assumed to be non-negative.
 181  */
 182 static int collect_fixed_sign_terms(__isl_take isl_term *term, void *user)
 183 {
 184         struct isl_fixed_sign_data *data = (struct isl_fixed_sign_data *)user;
 185         isl_int n;
 186         int i;
 187         int sign;
 188         unsigned nparam;
 189         unsigned nvar;
 190
 191         if (!term)
 192                 return -1;
 193
 194         nparam = isl_term_dim(term, isl_dim_param);
 195         nvar = isl_term_dim(term, isl_dim_set);
 196
 197         isl_int_init(n);
 198
 199         isl_term_get_num(term, &n);
 200
 201         sign = isl_int_sgn(n);
 202         for (i = 0; i < nparam; ++i) {
 203                 if (data->signs[i] > 0)
 204                         continue;
 205                 if (isl_term_get_exp(term, isl_dim_param, i) % 2)
 206                         sign = -sign;
 207         }
 208         for (i = 0; i < nvar; ++i) {
 209                 if (data->signs[nparam + i] > 0)
 210                         continue;
 211                 if (isl_term_get_exp(term, isl_dim_set, i) % 2)
 212                         sign = -sign;
 213         }
 214
 215         if (sign == data->sign) {
 216                 isl_qpolynomial *t = isl_qpolynomial_from_term(term);
 217
 218                 data->poly = isl_qpolynomial_add(data->poly, t);
 219         } else
 220                 isl_term_free(term);
 221
 222         isl_int_clear(n);
 223
 224         return 0;
 225 }
 226
 227 /* Construct and return a polynomial that consists of the terms
 228  * in "poly" that have sign "sign".  The integer divisions, if
 229  * any, are assumed to be non-negative.
 230  */
 231 __isl_give isl_qpolynomial *isl_qpolynomial_terms_of_sign(
 232         __isl_keep isl_qpolynomial *poly, int *signs, int sign)
 233 {
 234         isl_space *space;
 235         struct isl_fixed_sign_data data = { signs, sign };
 236
 237         space = isl_qpolynomial_get_domain_space(poly);
 238         data.poly = isl_qpolynomial_zero_on_domain(space);
 239
 240         if (isl_qpolynomial_foreach_term(poly, collect_fixed_sign_terms, &data) < 0)
 241                 goto error;
 242
 243         return data.poly;
 244 error:
 245         isl_qpolynomial_free(data.poly);
 246         return NULL;
 247 }
 248
 249 /* Helper function to add a guarded polynomial to either pwf_tight or pwf,
 250  * depending on whether the result has been determined to be tight.
 251  */
 252 static int add_guarded_poly(__isl_take isl_basic_set *bset,
 253         __isl_take isl_qpolynomial *poly, struct range_data *data)
 254 {
 255         enum isl_fold type = data->sign < 0 ? isl_fold_min : isl_fold_max;
 256         isl_set *set;
 257         isl_qpolynomial_fold *fold;
 258         isl_pw_qpolynomial_fold *pwf;
 259
 260         bset = isl_basic_set_params(bset);
 261         poly = isl_qpolynomial_project_domain_on_params(poly);
 262
 263         fold = isl_qpolynomial_fold_alloc(type, poly);
 264         set = isl_set_from_basic_set(bset);
 265         pwf = isl_pw_qpolynomial_fold_alloc(type, set, fold);
 266         if (data->tight)
 267                 data->pwf_tight = isl_pw_qpolynomial_fold_fold(
 268                                                 data->pwf_tight, pwf);
 269         else
 270                 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);
 271
 272         return 0;
 273 }
 274
 275 /* Given a lower and upper bound on the final variable and constraints
 276  * on the remaining variables where these bounds are active,
 277  * eliminate the variable from data->poly based on these bounds.
 278  * If the polynomial has been determined to be monotonic
 279  * in the variable, then simply plug in the appropriate bound.
 280  * If the current polynomial is tight and if this bound is integer,
 281  * then the result is still tight.  In all other cases, the results
 282  * may not be tight.
 283  * Otherwise, plug in the largest bound (in absolute value) in
 284  * the positive terms (if an upper bound is wanted) or the negative terms
 285  * (if a lower bounded is wanted) and the other bound in the other terms.
 286  *
 287  * If all variables have been eliminated, then record the result.
 288  * Ohterwise, recurse on the next variable.
 289  */
 290 static int propagate_on_bound_pair(__isl_take isl_constraint *lower,
 291         __isl_take isl_constraint *upper, __isl_take isl_basic_set *bset,
 292         void *user)
 293 {
 294         struct range_data *data = (struct range_data *)user;
 295         int save_tight = data->tight;
 296         isl_qpolynomial *poly;
 297         int r;
 298         unsigned nvar;
 299
 300         nvar = isl_basic_set_dim(bset, isl_dim_set);
 301
 302         if (data->monotonicity) {
 303                 isl_qpolynomial *sub;
 304                 isl_space *dim = isl_qpolynomial_get_domain_space(data->poly);
 305                 if (data->monotonicity * data->sign > 0) {
 306                         if (data->tight)
 307                                 data->tight = bound_is_integer(upper, nvar);
 308                         sub = bound2poly(upper, dim, nvar, 1);
 309                         isl_constraint_free(lower);
 310                 } else {
 311                         if (data->tight)
 312                                 data->tight = bound_is_integer(lower, nvar);
 313                         sub = bound2poly(lower, dim, nvar, -1);
 314                         isl_constraint_free(upper);
 315                 }
 316                 poly = isl_qpolynomial_copy(data->poly);
 317                 poly = isl_qpolynomial_substitute(poly, isl_dim_in, nvar, 1, &sub);
 318                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, nvar, 1);
 319
 320                 isl_qpolynomial_free(sub);
 321         } else {
 322                 isl_qpolynomial *l, *u;
 323                 isl_qpolynomial *pos, *neg;
 324                 isl_space *dim = isl_qpolynomial_get_domain_space(data->poly);
 325                 unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 326                 int sign = data->sign * data->signs[nparam + nvar];
 327
 328                 data->tight = 0;
 329
 330                 u = bound2poly(upper, isl_space_copy(dim), nvar, 1);
 331                 l = bound2poly(lower, dim, nvar, -1);
 332
 333                 pos = isl_qpolynomial_terms_of_sign(data->poly, data->signs, sign);
 334                 neg = isl_qpolynomial_terms_of_sign(data->poly, data->signs, -sign);
 335
 336                 pos = isl_qpolynomial_substitute(pos, isl_dim_in, nvar, 1, &u);
 337                 neg = isl_qpolynomial_substitute(neg, isl_dim_in, nvar, 1, &l);
 338
 339                 poly = isl_qpolynomial_add(pos, neg);
 340                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, nvar, 1);
 341
 342                 isl_qpolynomial_free(u);
 343                 isl_qpolynomial_free(l);
 344         }
 345
 346         if (isl_basic_set_dim(bset, isl_dim_set) == 0)
 347                 r = add_guarded_poly(bset, poly, data);
 348         else
 349                 r = propagate_on_domain(bset, poly, data);
 350
 351         data->tight = save_tight;
 352
 353         return r;
 354 }
 355
 356 /* Recursively perform range propagation on the polynomial "poly"
 357  * defined over the basic set "bset" and collect the results in "data".
 358  */
 359 static int propagate_on_domain(__isl_take isl_basic_set *bset,
 360         __isl_take isl_qpolynomial *poly, struct range_data *data)
 361 {
 362         isl_ctx *ctx;
 363         isl_qpolynomial *save_poly = data->poly;
 364         int save_monotonicity = data->monotonicity;
 365         unsigned d;
 366
 367         if (!bset || !poly)
 368                 goto error;
 369
 370         ctx = isl_basic_set_get_ctx(bset);
 371         d = isl_basic_set_dim(bset, isl_dim_set);
 372         isl_assert(ctx, d >= 1, goto error);
 373
 374         if (isl_qpolynomial_is_cst(poly, NULL, NULL)) {
 375                 bset = isl_basic_set_project_out(bset, isl_dim_set, 0, d);
 376                 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, d);
 377                 return add_guarded_poly(bset, poly, data);
 378         }
 379
 380         if (data->test_monotonicity)
 381                 data->monotonicity = monotonicity(bset, poly, data);
 382         else
 383                 data->monotonicity = 0;
 384         if (data->monotonicity < -1)
 385                 goto error;
 386
 387         data->poly = poly;
 388         if (isl_basic_set_foreach_bound_pair(bset, isl_dim_set, d - 1,
 389                                             &propagate_on_bound_pair, data) < 0)
 390                 goto error;
 391
 392         isl_basic_set_free(bset);
 393         isl_qpolynomial_free(poly);
 394         data->monotonicity = save_monotonicity;
 395         data->poly = save_poly;
 396
 397         return 0;
 398 error:
 399         isl_basic_set_free(bset);
 400         isl_qpolynomial_free(poly);
 401         data->monotonicity = save_monotonicity;
 402         data->poly = save_poly;
 403         return -1;
 404 }
 405
 406 static int basic_guarded_poly_bound(__isl_take isl_basic_set *bset, void *user)
 407 {
 408         struct range_data *data = (struct range_data *)user;
 409         isl_ctx *ctx;
 410         unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 411         unsigned dim = isl_basic_set_dim(bset, isl_dim_set);
 412         int r;
 413
 414         data->signs = NULL;
 415
 416         ctx = isl_basic_set_get_ctx(bset);
 417         data->signs = isl_alloc_array(ctx, int,
 418                                         isl_basic_set_dim(bset, isl_dim_all));
 419
 420         if (isl_basic_set_dims_get_sign(bset, isl_dim_set, 0, dim,
 421                                         data->signs + nparam) < 0)
 422                 goto error;
 423         if (isl_basic_set_dims_get_sign(bset, isl_dim_param, 0, nparam,
 424                                         data->signs) < 0)
 425                 goto error;
 426
 427         r = propagate_on_domain(bset, isl_qpolynomial_copy(data->poly), data);
 428
 429         free(data->signs);
 430
 431         return r;
 432 error:
 433         free(data->signs);
 434         isl_basic_set_free(bset);
 435         return -1;
 436 }
 437
 438 static int qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
 439         __isl_take isl_qpolynomial *poly, struct range_data *data)
 440 {
 441         unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
 442         unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
 443         isl_set *set;
 444
 445         if (!bset)
 446                 goto error;
 447
 448         if (nvar == 0)
 449                 return add_guarded_poly(bset, poly, data);
 450
 451         set = isl_set_from_basic_set(bset);
 452         set = isl_set_split_dims(set, isl_dim_param, 0, nparam);
 453         set = isl_set_split_dims(set, isl_dim_set, 0, nvar);
 454
 455         data->poly = poly;
 456
 457         data->test_monotonicity = 1;
 458         if (isl_set_foreach_basic_set(set, &basic_guarded_poly_bound, data) < 0)
 459                 goto error;
 460
 461         isl_set_free(set);
 462         isl_qpolynomial_free(poly);
 463
 464         return 0;
 465 error:
 466         isl_set_free(set);
 467         isl_qpolynomial_free(poly);
 468         return -1;
 469 }
 470
 471 int isl_qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
 472         __isl_take isl_qpolynomial *poly, struct isl_bound *bound)
 473 {
 474         struct range_data data;
 475         int r;
 476
 477         data.pwf = bound->pwf;
 478         data.pwf_tight = bound->pwf_tight;
 479         data.tight = bound->check_tight;
 480         if (bound->type == isl_fold_min)
 481                 data.sign = -1;
 482         else
 483                 data.sign = 1;
 484
 485         r = qpolynomial_bound_on_domain_range(bset, poly, &data);
 486
 487         bound->pwf = data.pwf;
 488         bound->pwf_tight = data.pwf_tight;
 489
 490         return r;
 491 }