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_bernstein.h> 26 27struct bernstein_data { 28 enum isl_fold type; 29 isl_qpolynomial *poly; 30 int check_tight; 31 32 isl_cell *cell; 33 34 isl_qpolynomial_fold *fold; 35 isl_qpolynomial_fold *fold_tight; 36 isl_pw_qpolynomial_fold *pwf; 37 isl_pw_qpolynomial_fold *pwf_tight; 38}; 39 40static int vertex_is_integral(__isl_keep isl_basic_set *vertex) 41{ 42 unsigned nvar; 43 unsigned nparam; 44 int i; 45 46 nvar = isl_basic_set_dim(vertex, isl_dim_set); 47 nparam = isl_basic_set_dim(vertex, isl_dim_param); 48 for (i = 0; i < nvar; ++i) { 49 int r = nvar - 1 - i; 50 if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) && 51 !isl_int_is_negone(vertex->eq[r][1 + nparam + i])) 52 return 0; 53 } 54 55 return 1; 56} 57 58static __isl_give isl_qpolynomial *vertex_coordinate( 59 __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *dim) 60{ 61 unsigned nvar; 62 unsigned nparam; 63 int r; 64 isl_int denom; 65 isl_qpolynomial *v; 66 67 nvar = isl_basic_set_dim(vertex, isl_dim_set); 68 nparam = isl_basic_set_dim(vertex, isl_dim_param); 69 r = nvar - 1 - i; 70 71 isl_int_init(denom); 72 isl_int_set(denom, vertex->eq[r][1 + nparam + i]); 73 isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error); 74 75 if (isl_int_is_pos(denom)) 76 isl_seq_neg(vertex->eq[r], vertex->eq[r], 77 1 + isl_basic_set_total_dim(vertex)); 78 else 79 isl_int_neg(denom, denom); 80 81 v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom); 82 isl_int_clear(denom); 83 84 return v; 85error: 86 isl_space_free(dim); 87 isl_int_clear(denom); 88 return NULL; 89} 90 91/* Check whether the bound associated to the selection "k" is tight, 92 * which is the case if we select exactly one vertex and if that vertex 93 * is integral for all values of the parameters. 94 */ 95static int is_tight(int *k, int n, int d, isl_cell *cell) 96{ 97 int i; 98 99 for (i = 0; i < n; ++i) { 100 int v; 101 if (k[i] != d) { 102 if (k[i]) 103 return 0; 104 continue; 105 } 106 v = cell->ids[n - 1 - i]; 107 return vertex_is_integral(cell->vertices->v[v].vertex); 108 } 109 110 return 0; 111} 112 113static void add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom, 114 int *k, int n, int d, struct bernstein_data *data) 115{ 116 isl_qpolynomial_fold *fold; 117 118 fold = isl_qpolynomial_fold_alloc(data->type, b); 119 120 if (data->check_tight && is_tight(k, n, d, data->cell)) 121 data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom, 122 data->fold_tight, fold); 123 else 124 data->fold = isl_qpolynomial_fold_fold_on_domain(dom, 125 data->fold, fold); 126} 127 128/* Extract the coefficients of the Bernstein base polynomials and store 129 * them in data->fold and data->fold_tight. 130 * 131 * In particular, the coefficient of each monomial 132 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding 133 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]! 134 * 135 * c[i] contains the coefficient of the selected powers of the first i+1 vars. 136 * multinom[i] contains the partial multinomial coefficient. 137 */ 138static void extract_coefficients(isl_qpolynomial *poly, 139 __isl_keep isl_set *dom, struct bernstein_data *data) 140{ 141 int i; 142 int d; 143 int n; 144 isl_ctx *ctx; 145 isl_qpolynomial **c = NULL; 146 int *k = NULL; 147 int *left = NULL; 148 isl_vec *multinom = NULL; 149 150 if (!poly) 151 return; 152 153 ctx = isl_qpolynomial_get_ctx(poly); 154 n = isl_qpolynomial_dim(poly, isl_dim_in); 155 d = isl_qpolynomial_degree(poly); 156 isl_assert(ctx, n >= 2, return); 157 158 c = isl_calloc_array(ctx, isl_qpolynomial *, n); 159 k = isl_alloc_array(ctx, int, n); 160 left = isl_alloc_array(ctx, int, n); 161 multinom = isl_vec_alloc(ctx, n); 162 if (!c || !k || !left || !multinom) 163 goto error; 164 165 isl_int_set_si(multinom->el[0], 1); 166 for (k[0] = d; k[0] >= 0; --k[0]) { 167 int i = 1; 168 isl_qpolynomial_free(c[0]); 169 c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]); 170 left[0] = d - k[0]; 171 k[1] = -1; 172 isl_int_set(multinom->el[1], multinom->el[0]); 173 while (i > 0) { 174 if (i == n - 1) { 175 int j; 176 isl_space *dim; 177 isl_qpolynomial *b; 178 isl_qpolynomial *f; 179 for (j = 2; j <= left[i - 1]; ++j) 180 isl_int_divexact_ui(multinom->el[i], 181 multinom->el[i], j); 182 b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in, 183 n - 1 - i, left[i - 1]); 184 b = isl_qpolynomial_project_domain_on_params(b); 185 dim = isl_qpolynomial_get_domain_space(b); 186 f = isl_qpolynomial_rat_cst_on_domain(dim, ctx->one, 187 multinom->el[i]); 188 b = isl_qpolynomial_mul(b, f); 189 k[n - 1] = left[n - 2]; 190 add_fold(b, dom, k, n, d, data); 191 --i; 192 continue; 193 } 194 if (k[i] >= left[i - 1]) { 195 --i; 196 continue; 197 } 198 ++k[i]; 199 if (k[i]) 200 isl_int_divexact_ui(multinom->el[i], 201 multinom->el[i], k[i]); 202 isl_qpolynomial_free(c[i]); 203 c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in, 204 n - 1 - i, k[i]); 205 left[i] = left[i - 1] - k[i]; 206 k[i + 1] = -1; 207 isl_int_set(multinom->el[i + 1], multinom->el[i]); 208 ++i; 209 } 210 isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]); 211 } 212 213 for (i = 0; i < n; ++i) 214 isl_qpolynomial_free(c[i]); 215 216 isl_vec_free(multinom); 217 free(left); 218 free(k); 219 free(c); 220 return; 221error: 222 isl_vec_free(multinom); 223 free(left); 224 free(k); 225 if (c) 226 for (i = 0; i < n; ++i) 227 isl_qpolynomial_free(c[i]); 228 free(c); 229 return; 230} 231 232/* Perform bernstein expansion on the parametric vertices that are active 233 * on "cell". 234 * 235 * data->poly has been homogenized in the calling function. 236 * 237 * We plug in the barycentric coordinates for the set variables 238 * 239 * \vec x = \sum_i \alpha_i v_i(\vec p) 240 * 241 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension. 242 * Next, we extract the coefficients of the Bernstein base polynomials. 243 */ 244static int bernstein_coefficients_cell(__isl_take isl_cell *cell, void *user) 245{ 246 int i, j; 247 struct bernstein_data *data = (struct bernstein_data *)user; 248 isl_space *dim_param; 249 isl_space *dim_dst; 250 isl_qpolynomial *poly = data->poly; 251 unsigned nvar; 252 int n_vertices; 253 isl_qpolynomial **subs; 254 isl_pw_qpolynomial_fold *pwf; 255 isl_set *dom; 256 isl_ctx *ctx; 257 258 if (!poly) 259 goto error; 260 261 nvar = isl_qpolynomial_dim(poly, isl_dim_in) - 1; 262 n_vertices = cell->n_vertices; 263 264 ctx = isl_qpolynomial_get_ctx(poly); 265 if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate) 266 return isl_cell_foreach_simplex(cell, 267 &bernstein_coefficients_cell, user); 268 269 subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar); 270 if (!subs) 271 goto error; 272 273 dim_param = isl_basic_set_get_space(cell->dom); 274 dim_dst = isl_qpolynomial_get_domain_space(poly); 275 dim_dst = isl_space_add_dims(dim_dst, isl_dim_set, n_vertices); 276 277 for (i = 0; i < 1 + nvar; ++i) 278 subs[i] = isl_qpolynomial_zero_on_domain(isl_space_copy(dim_dst)); 279 280 for (i = 0; i < n_vertices; ++i) { 281 isl_qpolynomial *c; 282 c = isl_qpolynomial_var_on_domain(isl_space_copy(dim_dst), isl_dim_set, 283 1 + nvar + i); 284 for (j = 0; j < nvar; ++j) { 285 int k = cell->ids[i]; 286 isl_qpolynomial *v; 287 v = vertex_coordinate(cell->vertices->v[k].vertex, j, 288 isl_space_copy(dim_param)); 289 v = isl_qpolynomial_add_dims(v, isl_dim_in, 290 1 + nvar + n_vertices); 291 v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c)); 292 subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v); 293 } 294 subs[0] = isl_qpolynomial_add(subs[0], c); 295 } 296 isl_space_free(dim_dst); 297 298 poly = isl_qpolynomial_copy(poly); 299 300 poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices); 301 poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs); 302 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar); 303 304 data->cell = cell; 305 dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom)); 306 data->fold = isl_qpolynomial_fold_empty(data->type, isl_space_copy(dim_param)); 307 data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param); 308 extract_coefficients(poly, dom, data); 309 310 pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom), 311 data->fold); 312 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf); 313 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight); 314 data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf); 315 316 isl_qpolynomial_free(poly); 317 isl_cell_free(cell); 318 for (i = 0; i < 1 + nvar; ++i) 319 isl_qpolynomial_free(subs[i]); 320 free(subs); 321 return 0; 322error: 323 isl_cell_free(cell); 324 return -1; 325} 326 327/* Base case of applying bernstein expansion. 328 * 329 * We compute the chamber decomposition of the parametric polytope "bset" 330 * and then perform bernstein expansion on the parametric vertices 331 * that are active on each chamber. 332 */ 333static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base( 334 __isl_take isl_basic_set *bset, 335 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) 336{ 337 unsigned nvar; 338 isl_space *dim; 339 isl_pw_qpolynomial_fold *pwf; 340 isl_vertices *vertices; 341 int covers; 342 343 nvar = isl_basic_set_dim(bset, isl_dim_set); 344 if (nvar == 0) { 345 isl_set *dom; 346 isl_qpolynomial_fold *fold; 347 348 fold = isl_qpolynomial_fold_alloc(data->type, poly); 349 dom = isl_set_from_basic_set(bset); 350 if (tight) 351 *tight = 1; 352 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold); 353 return isl_pw_qpolynomial_fold_project_domain_on_params(pwf); 354 } 355 356 if (isl_qpolynomial_is_zero(poly)) { 357 isl_set *dom; 358 isl_qpolynomial_fold *fold; 359 fold = isl_qpolynomial_fold_alloc(data->type, poly); 360 dom = isl_set_from_basic_set(bset); 361 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold); 362 if (tight) 363 *tight = 1; 364 return isl_pw_qpolynomial_fold_project_domain_on_params(pwf); 365 } 366 367 dim = isl_basic_set_get_space(bset); 368 dim = isl_space_params(dim); 369 dim = isl_space_from_domain(dim); 370 dim = isl_space_add_dims(dim, isl_dim_set, 1); 371 data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(dim), data->type); 372 data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type); 373 data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly)); 374 vertices = isl_basic_set_compute_vertices(bset); 375 isl_vertices_foreach_disjoint_cell(vertices, 376 &bernstein_coefficients_cell, data); 377 isl_vertices_free(vertices); 378 isl_qpolynomial_free(data->poly); 379 380 isl_basic_set_free(bset); 381 isl_qpolynomial_free(poly); 382 383 covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf); 384 if (covers < 0) 385 goto error; 386 387 if (tight) 388 *tight = covers; 389 390 if (covers) { 391 isl_pw_qpolynomial_fold_free(data->pwf); 392 return data->pwf_tight; 393 } 394 395 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight); 396 397 return data->pwf; 398error: 399 isl_pw_qpolynomial_fold_free(data->pwf_tight); 400 isl_pw_qpolynomial_fold_free(data->pwf); 401 return NULL; 402} 403 404/* Apply bernstein expansion recursively by working in on len[i] 405 * set variables at a time, with i ranging from n_group - 1 to 0. 406 */ 407static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive( 408 __isl_take isl_pw_qpolynomial *pwqp, 409 int n_group, int *len, struct bernstein_data *data, int *tight) 410{ 411 int i; 412 unsigned nparam; 413 unsigned nvar; 414 isl_pw_qpolynomial_fold *pwf; 415 416 if (!pwqp) 417 return NULL; 418 419 nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param); 420 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in); 421 422 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam, 423 isl_dim_in, 0, nvar - len[n_group - 1]); 424 pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight); 425 426 for (i = n_group - 2; i >= 0; --i) { 427 nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param); 428 pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0, 429 isl_dim_param, nparam - len[i], len[i]); 430 if (tight && !*tight) 431 tight = NULL; 432 pwf = isl_pw_qpolynomial_fold_bound(pwf, tight); 433 } 434 435 return pwf; 436} 437 438static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors( 439 __isl_take isl_basic_set *bset, 440 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) 441{ 442 isl_factorizer *f; 443 isl_set *set; 444 isl_pw_qpolynomial *pwqp; 445 isl_pw_qpolynomial_fold *pwf; 446 447 f = isl_basic_set_factorizer(bset); 448 if (!f) 449 goto error; 450 if (f->n_group == 0) { 451 isl_factorizer_free(f); 452 return bernstein_coefficients_base(bset, poly, data, tight); 453 } 454 455 set = isl_set_from_basic_set(bset); 456 pwqp = isl_pw_qpolynomial_alloc(set, poly); 457 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph)); 458 459 pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data, 460 tight); 461 462 isl_factorizer_free(f); 463 464 return pwf; 465error: 466 isl_basic_set_free(bset); 467 isl_qpolynomial_free(poly); 468 return NULL; 469} 470 471static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive( 472 __isl_take isl_basic_set *bset, 473 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) 474{ 475 int i; 476 int *len; 477 unsigned nvar; 478 isl_pw_qpolynomial_fold *pwf; 479 isl_set *set; 480 isl_pw_qpolynomial *pwqp; 481 482 if (!bset || !poly) 483 goto error; 484 485 nvar = isl_basic_set_dim(bset, isl_dim_set); 486 487 len = isl_alloc_array(bset->ctx, int, nvar); 488 if (nvar && !len) 489 goto error; 490 491 for (i = 0; i < nvar; ++i) 492 len[i] = 1; 493 494 set = isl_set_from_basic_set(bset); 495 pwqp = isl_pw_qpolynomial_alloc(set, poly); 496 497 pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight); 498 499 free(len); 500 501 return pwf; 502error: 503 isl_basic_set_free(bset); 504 isl_qpolynomial_free(poly); 505 return NULL; 506} 507 508/* Compute a bound on the polynomial defined over the parametric polytope 509 * using bernstein expansion and store the result 510 * in bound->pwf and bound->pwf_tight. 511 * 512 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if 513 * the polytope can be factorized and apply bernstein expansion recursively 514 * on the factors. 515 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply 516 * bernstein expansion recursively on each dimension. 517 * Otherwise, we apply bernstein expansion on the entire polytope. 518 */ 519int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set *bset, 520 __isl_take isl_qpolynomial *poly, struct isl_bound *bound) 521{ 522 struct bernstein_data data; 523 isl_pw_qpolynomial_fold *pwf; 524 unsigned nvar; 525 int tight = 0; 526 int *tp = bound->check_tight ? &tight : NULL; 527 528 if (!bset || !poly) 529 goto error; 530 531 data.type = bound->type; 532 data.check_tight = bound->check_tight; 533 534 nvar = isl_basic_set_dim(bset, isl_dim_set); 535 536 if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS) 537 pwf = bernstein_coefficients_factors(bset, poly, &data, tp); 538 else if (nvar > 1 && 539 (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS)) 540 pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp); 541 else 542 pwf = bernstein_coefficients_base(bset, poly, &data, tp); 543 544 if (tight) 545 bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf); 546 else 547 bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf); 548 549 return 0; 550error: 551 isl_basic_set_free(bset); 552 isl_qpolynomial_free(poly); 553 return -1; 554} 555