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