584 lines · c
1/*2 * Copyright 2006-2007 Universiteit Leiden3 * Copyright 2008-2009 Katholieke Universiteit Leuven4 * Copyright 2010 INRIA Saclay5 *6 * Use of this software is governed by the MIT license7 *8 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science,9 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands10 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A,11 * B-3001 Leuven, Belgium12 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,13 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France14 */15 16#include <isl_ctx_private.h>17#include <isl_map_private.h>18#include <isl/set.h>19#include <isl_seq.h>20#include <isl_morph.h>21#include <isl_factorization.h>22#include <isl_vertices_private.h>23#include <isl_polynomial_private.h>24#include <isl_options_private.h>25#include <isl_vec_private.h>26#include <isl_bernstein.h>27 28struct bernstein_data {29 enum isl_fold type;30 isl_qpolynomial *poly;31 int check_tight;32 33 isl_cell *cell;34 35 isl_qpolynomial_fold *fold;36 isl_qpolynomial_fold *fold_tight;37 isl_pw_qpolynomial_fold *pwf;38 isl_pw_qpolynomial_fold *pwf_tight;39};40 41static isl_bool vertex_is_integral(__isl_keep isl_basic_set *vertex)42{43 isl_size nvar;44 isl_size nparam;45 int i;46 47 nvar = isl_basic_set_dim(vertex, isl_dim_set);48 nparam = isl_basic_set_dim(vertex, isl_dim_param);49 if (nvar < 0 || nparam < 0)50 return isl_bool_error;51 for (i = 0; i < nvar; ++i) {52 int r = nvar - 1 - i;53 if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&54 !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))55 return isl_bool_false;56 }57 58 return isl_bool_true;59}60 61static __isl_give isl_qpolynomial *vertex_coordinate(62 __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *space)63{64 isl_size nvar;65 isl_size nparam;66 isl_size total;67 int r;68 isl_int denom;69 isl_qpolynomial *v;70 71 isl_int_init(denom);72 73 nvar = isl_basic_set_dim(vertex, isl_dim_set);74 nparam = isl_basic_set_dim(vertex, isl_dim_param);75 total = isl_basic_set_dim(vertex, isl_dim_all);76 if (nvar < 0 || nparam < 0 || total < 0)77 goto error;78 r = nvar - 1 - i;79 80 isl_int_set(denom, vertex->eq[r][1 + nparam + i]);81 isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);82 83 if (isl_int_is_pos(denom))84 isl_seq_neg(vertex->eq[r], vertex->eq[r], 1 + total);85 else86 isl_int_neg(denom, denom);87 88 v = isl_qpolynomial_from_affine(space, vertex->eq[r], denom);89 isl_int_clear(denom);90 91 return v;92error:93 isl_space_free(space);94 isl_int_clear(denom);95 return NULL;96}97 98/* Check whether the bound associated to the selection "k" is tight,99 * which is the case if we select exactly one vertex (i.e., one of the100 * exponents in "k" is exactly "d") and if that vertex101 * is integral for all values of the parameters.102 *103 * If the degree "d" is zero, then there are no exponents.104 * Since the polynomial is a constant expression in this case,105 * the bound is necessarily tight.106 */107static isl_bool is_tight(int *k, int n, int d, isl_cell *cell)108{109 int i;110 111 if (d == 0)112 return isl_bool_true;113 114 for (i = 0; i < n; ++i) {115 int v;116 if (!k[i])117 continue;118 if (k[i] != d)119 return isl_bool_false;120 v = cell->ids[n - 1 - i];121 return vertex_is_integral(cell->vertices->v[v].vertex);122 }123 124 return isl_bool_false;125}126 127static isl_stat add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom,128 int *k, int n, int d, struct bernstein_data *data)129{130 isl_qpolynomial_fold *fold;131 isl_bool tight;132 133 fold = isl_qpolynomial_fold_alloc(data->type, b);134 135 tight = isl_bool_false;136 if (data->check_tight)137 tight = is_tight(k, n, d, data->cell);138 if (tight < 0)139 return isl_stat_error;140 if (tight)141 data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom,142 data->fold_tight, fold);143 else144 data->fold = isl_qpolynomial_fold_fold_on_domain(dom,145 data->fold, fold);146 return isl_stat_ok;147}148 149/* Extract the coefficients of the Bernstein base polynomials and store150 * them in data->fold and data->fold_tight.151 *152 * In particular, the coefficient of each monomial153 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding154 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]!155 *156 * c[i] contains the coefficient of the selected powers of the first i+1 vars.157 * multinom[i] contains the partial multinomial coefficient.158 */159static isl_stat extract_coefficients(isl_qpolynomial *poly,160 __isl_keep isl_set *dom, struct bernstein_data *data)161{162 int i;163 int d;164 isl_size n;165 isl_ctx *ctx;166 isl_qpolynomial **c = NULL;167 int *k = NULL;168 int *left = NULL;169 isl_vec *multinom = NULL;170 171 n = isl_qpolynomial_dim(poly, isl_dim_in);172 if (n < 0)173 return isl_stat_error;174 175 ctx = isl_qpolynomial_get_ctx(poly);176 d = isl_qpolynomial_degree(poly);177 isl_assert(ctx, n >= 2, return isl_stat_error);178 179 c = isl_calloc_array(ctx, isl_qpolynomial *, n);180 k = isl_alloc_array(ctx, int, n);181 left = isl_alloc_array(ctx, int, n);182 multinom = isl_vec_alloc(ctx, n);183 if (!c || !k || !left || !multinom)184 goto error;185 186 isl_int_set_si(multinom->el[0], 1);187 for (k[0] = d; k[0] >= 0; --k[0]) {188 int i = 1;189 isl_qpolynomial_free(c[0]);190 c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]);191 left[0] = d - k[0];192 k[1] = -1;193 isl_int_set(multinom->el[1], multinom->el[0]);194 while (i > 0) {195 if (i == n - 1) {196 int j;197 isl_space *space;198 isl_qpolynomial *b;199 isl_qpolynomial *f;200 for (j = 2; j <= left[i - 1]; ++j)201 isl_int_divexact_ui(multinom->el[i],202 multinom->el[i], j);203 b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,204 n - 1 - i, left[i - 1]);205 b = isl_qpolynomial_project_domain_on_params(b);206 space = isl_qpolynomial_get_domain_space(b);207 f = isl_qpolynomial_rat_cst_on_domain(space,208 ctx->one, multinom->el[i]);209 b = isl_qpolynomial_mul(b, f);210 k[n - 1] = left[n - 2];211 if (add_fold(b, dom, k, n, d, data) < 0)212 goto error;213 --i;214 continue;215 }216 if (k[i] >= left[i - 1]) {217 --i;218 continue;219 }220 ++k[i];221 if (k[i])222 isl_int_divexact_ui(multinom->el[i],223 multinom->el[i], k[i]);224 isl_qpolynomial_free(c[i]);225 c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,226 n - 1 - i, k[i]);227 left[i] = left[i - 1] - k[i];228 k[i + 1] = -1;229 isl_int_set(multinom->el[i + 1], multinom->el[i]);230 ++i;231 }232 isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);233 }234 235 for (i = 0; i < n; ++i)236 isl_qpolynomial_free(c[i]);237 238 isl_vec_free(multinom);239 free(left);240 free(k);241 free(c);242 return isl_stat_ok;243error:244 isl_vec_free(multinom);245 free(left);246 free(k);247 if (c)248 for (i = 0; i < n; ++i)249 isl_qpolynomial_free(c[i]);250 free(c);251 return isl_stat_error;252}253 254/* Perform bernstein expansion on the parametric vertices that are active255 * on "cell".256 *257 * data->poly has been homogenized in the calling function.258 *259 * We plug in the barycentric coordinates for the set variables260 *261 * \vec x = \sum_i \alpha_i v_i(\vec p)262 *263 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.264 * Next, we extract the coefficients of the Bernstein base polynomials.265 */266static isl_stat bernstein_coefficients_cell(__isl_take isl_cell *cell,267 void *user)268{269 int i, j;270 struct bernstein_data *data = (struct bernstein_data *)user;271 isl_space *space_param;272 isl_space *space_dst;273 isl_qpolynomial *poly = data->poly;274 isl_size n_in;275 unsigned nvar;276 int n_vertices;277 isl_qpolynomial **subs;278 isl_pw_qpolynomial_fold *pwf;279 isl_set *dom;280 isl_ctx *ctx;281 282 n_in = isl_qpolynomial_dim(poly, isl_dim_in);283 if (n_in < 0)284 goto error;285 286 nvar = n_in - 1;287 n_vertices = cell->n_vertices;288 289 ctx = isl_qpolynomial_get_ctx(poly);290 if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate)291 return isl_cell_foreach_simplex(cell,292 &bernstein_coefficients_cell, user);293 294 subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar);295 if (!subs)296 goto error;297 298 space_param = isl_basic_set_get_space(cell->dom);299 space_dst = isl_qpolynomial_get_domain_space(poly);300 space_dst = isl_space_add_dims(space_dst, isl_dim_set, n_vertices);301 302 for (i = 0; i < 1 + nvar; ++i)303 subs[i] =304 isl_qpolynomial_zero_on_domain(isl_space_copy(space_dst));305 306 for (i = 0; i < n_vertices; ++i) {307 isl_qpolynomial *c;308 c = isl_qpolynomial_var_on_domain(isl_space_copy(space_dst),309 isl_dim_set, 1 + nvar + i);310 for (j = 0; j < nvar; ++j) {311 int k = cell->ids[i];312 isl_qpolynomial *v;313 v = vertex_coordinate(cell->vertices->v[k].vertex, j,314 isl_space_copy(space_param));315 v = isl_qpolynomial_add_dims(v, isl_dim_in,316 1 + nvar + n_vertices);317 v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c));318 subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v);319 }320 subs[0] = isl_qpolynomial_add(subs[0], c);321 }322 isl_space_free(space_dst);323 324 poly = isl_qpolynomial_copy(poly);325 326 poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices);327 poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs);328 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar);329 330 data->cell = cell;331 dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom));332 data->fold = isl_qpolynomial_fold_empty(data->type,333 isl_space_copy(space_param));334 data->fold_tight = isl_qpolynomial_fold_empty(data->type, space_param);335 if (extract_coefficients(poly, dom, data) < 0) {336 data->fold = isl_qpolynomial_fold_free(data->fold);337 data->fold_tight = isl_qpolynomial_fold_free(data->fold_tight);338 }339 340 pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),341 data->fold);342 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);343 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight);344 data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf);345 346 isl_qpolynomial_free(poly);347 isl_cell_free(cell);348 for (i = 0; i < 1 + nvar; ++i)349 isl_qpolynomial_free(subs[i]);350 free(subs);351 return isl_stat_ok;352error:353 isl_cell_free(cell);354 return isl_stat_error;355}356 357/* Base case of applying bernstein expansion.358 *359 * We compute the chamber decomposition of the parametric polytope "bset"360 * and then perform bernstein expansion on the parametric vertices361 * that are active on each chamber.362 *363 * If the polynomial does not depend on the set variables364 * (and in particular if the number of set variables is zero)365 * then the bound is equal to the polynomial and366 * no actual bernstein expansion needs to be performed.367 */368static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base(369 __isl_take isl_basic_set *bset,370 __isl_take isl_qpolynomial *poly, struct bernstein_data *data,371 isl_bool *tight)372{373 int degree;374 isl_size nvar;375 isl_space *space;376 isl_vertices *vertices;377 isl_bool covers;378 379 nvar = isl_basic_set_dim(bset, isl_dim_set);380 if (nvar < 0)381 bset = isl_basic_set_free(bset);382 if (nvar == 0)383 return isl_qpolynomial_cst_bound(bset, poly, data->type, tight);384 385 degree = isl_qpolynomial_degree(poly);386 if (degree < -1)387 bset = isl_basic_set_free(bset);388 if (degree <= 0)389 return isl_qpolynomial_cst_bound(bset, poly, data->type, tight);390 391 space = isl_basic_set_get_space(bset);392 space = isl_space_params(space);393 space = isl_space_from_domain(space);394 space = isl_space_add_dims(space, isl_dim_set, 1);395 data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(space),396 data->type);397 data->pwf_tight = isl_pw_qpolynomial_fold_zero(space, data->type);398 data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly));399 vertices = isl_basic_set_compute_vertices(bset);400 if (isl_vertices_foreach_disjoint_cell(vertices,401 &bernstein_coefficients_cell, data) < 0)402 data->pwf = isl_pw_qpolynomial_fold_free(data->pwf);403 isl_vertices_free(vertices);404 isl_qpolynomial_free(data->poly);405 406 isl_basic_set_free(bset);407 isl_qpolynomial_free(poly);408 409 covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);410 if (covers < 0)411 goto error;412 413 if (tight)414 *tight = covers;415 416 if (covers) {417 isl_pw_qpolynomial_fold_free(data->pwf);418 return data->pwf_tight;419 }420 421 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);422 423 return data->pwf;424error:425 isl_pw_qpolynomial_fold_free(data->pwf_tight);426 isl_pw_qpolynomial_fold_free(data->pwf);427 return NULL;428}429 430/* Apply bernstein expansion recursively by working in on len[i]431 * set variables at a time, with i ranging from n_group - 1 to 0.432 */433static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive(434 __isl_take isl_pw_qpolynomial *pwqp,435 int n_group, int *len, struct bernstein_data *data, isl_bool *tight)436{437 int i;438 isl_size nparam;439 isl_size nvar;440 isl_pw_qpolynomial_fold *pwf;441 442 nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);443 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in);444 if (nparam < 0 || nvar < 0)445 goto error;446 447 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,448 isl_dim_in, 0, nvar - len[n_group - 1]);449 pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight);450 451 for (i = n_group - 2; i >= 0; --i) {452 nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);453 if (nparam < 0)454 return isl_pw_qpolynomial_fold_free(pwf);455 pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0,456 isl_dim_param, nparam - len[i], len[i]);457 if (tight && !*tight)458 tight = NULL;459 pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);460 }461 462 return pwf;463error:464 isl_pw_qpolynomial_free(pwqp);465 return NULL;466}467 468static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors(469 __isl_take isl_basic_set *bset,470 __isl_take isl_qpolynomial *poly, struct bernstein_data *data,471 isl_bool *tight)472{473 isl_factorizer *f;474 isl_set *set;475 isl_pw_qpolynomial *pwqp;476 isl_pw_qpolynomial_fold *pwf;477 478 f = isl_basic_set_factorizer(bset);479 if (!f)480 goto error;481 if (f->n_group == 0) {482 isl_factorizer_free(f);483 return bernstein_coefficients_base(bset, poly, data, tight);484 }485 486 set = isl_set_from_basic_set(bset);487 pwqp = isl_pw_qpolynomial_alloc(set, poly);488 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph));489 490 pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,491 tight);492 493 isl_factorizer_free(f);494 495 return pwf;496error:497 isl_basic_set_free(bset);498 isl_qpolynomial_free(poly);499 return NULL;500}501 502static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive(503 __isl_take isl_basic_set *bset,504 __isl_take isl_qpolynomial *poly, struct bernstein_data *data,505 isl_bool *tight)506{507 int i;508 int *len;509 isl_size nvar;510 isl_pw_qpolynomial_fold *pwf;511 isl_set *set;512 isl_pw_qpolynomial *pwqp;513 514 nvar = isl_basic_set_dim(bset, isl_dim_set);515 if (nvar < 0 || !poly)516 goto error;517 518 len = isl_alloc_array(bset->ctx, int, nvar);519 if (nvar && !len)520 goto error;521 522 for (i = 0; i < nvar; ++i)523 len[i] = 1;524 525 set = isl_set_from_basic_set(bset);526 pwqp = isl_pw_qpolynomial_alloc(set, poly);527 528 pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);529 530 free(len);531 532 return pwf;533error:534 isl_basic_set_free(bset);535 isl_qpolynomial_free(poly);536 return NULL;537}538 539/* Compute a bound on the polynomial defined over the parametric polytope540 * using bernstein expansion and store the result541 * in bound->pwf and bound->pwf_tight.542 *543 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if544 * the polytope can be factorized and apply bernstein expansion recursively545 * on the factors.546 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply547 * bernstein expansion recursively on each dimension.548 * Otherwise, we apply bernstein expansion on the entire polytope.549 */550isl_stat isl_qpolynomial_bound_on_domain_bernstein(551 __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *poly,552 struct isl_bound *bound)553{554 struct bernstein_data data;555 isl_pw_qpolynomial_fold *pwf;556 isl_size nvar;557 isl_bool tight = isl_bool_false;558 isl_bool *tp = bound->check_tight ? &tight : NULL;559 560 nvar = isl_basic_set_dim(bset, isl_dim_set);561 if (nvar < 0 || !poly)562 goto error;563 564 data.type = bound->type;565 data.check_tight = bound->check_tight;566 567 if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)568 pwf = bernstein_coefficients_factors(bset, poly, &data, tp);569 else if (nvar > 1 &&570 (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))571 pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);572 else573 pwf = bernstein_coefficients_base(bset, poly, &data, tp);574 575 if (tight)576 return isl_bound_add_tight(bound, pwf);577 else578 return isl_bound_add(bound, pwf);579error:580 isl_basic_set_free(bset);581 isl_qpolynomial_free(poly);582 return isl_stat_error;583}584