1/* 2 * Copyright 2006-2007 Universiteit Leiden 3 * Copyright 2008-2009 Katholieke Universiteit Leuven 4 * Copyright 2010 INRIA Saclay 5 * 6 * Use of this software is governed by the MIT license 7 * 8 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science, 9 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands 10 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A, 11 * B-3001 Leuven, Belgium 12 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, 13 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 14 */ 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 else 86 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 the 100 * exponents in "k" is exactly "d") and if that vertex 101 * 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 else 144 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 store 150 * them in data->fold and data->fold_tight. 151 * 152 * In particular, the coefficient of each monomial 153 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding 154 * 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 active 255 * on "cell". 256 * 257 * data->poly has been homogenized in the calling function. 258 * 259 * We plug in the barycentric coordinates for the set variables 260 * 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 vertices 361 * that are active on each chamber. 362 * 363 * If the polynomial does not depend on the set variables 364 * (and in particular if the number of set variables is zero) 365 * then the bound is equal to the polynomial and 366 * 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 polytope 540 * using bernstein expansion and store the result 541 * in bound->pwf and bound->pwf_tight. 542 * 543 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if 544 * the polytope can be factorized and apply bernstein expansion recursively 545 * on the factors. 546 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply 547 * 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 else 573 pwf = bernstein_coefficients_base(bset, poly, &data, tp); 574 575 if (tight) 576 return isl_bound_add_tight(bound, pwf); 577 else 578 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