diff options
Diffstat (limited to 'cloog-0.17.0/isl/isl_bernstein.c')
| -rw-r--r-- | cloog-0.17.0/isl/isl_bernstein.c | 554 |
1 files changed, 0 insertions, 554 deletions
diff --git a/cloog-0.17.0/isl/isl_bernstein.c b/cloog-0.17.0/isl/isl_bernstein.c deleted file mode 100644 index 0e528e5..0000000 --- a/cloog-0.17.0/isl/isl_bernstein.c +++ /dev/null @@ -1,554 +0,0 @@ -/* - * Copyright 2006-2007 Universiteit Leiden - * Copyright 2008-2009 Katholieke Universiteit Leuven - * Copyright 2010 INRIA Saclay - * - * Use of this software is governed by the GNU LGPLv2.1 license - * - * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science, - * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands - * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A, - * B-3001 Leuven, Belgium - * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, - * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France - */ - -#include <isl_ctx_private.h> -#include <isl_map_private.h> -#include <isl/set.h> -#include <isl/seq.h> -#include <isl_morph.h> -#include <isl_factorization.h> -#include <isl_vertices_private.h> -#include <isl_polynomial_private.h> -#include <isl_options_private.h> -#include <isl_bernstein.h> - -struct bernstein_data { - enum isl_fold type; - isl_qpolynomial *poly; - int check_tight; - - isl_cell *cell; - - isl_qpolynomial_fold *fold; - isl_qpolynomial_fold *fold_tight; - isl_pw_qpolynomial_fold *pwf; - isl_pw_qpolynomial_fold *pwf_tight; -}; - -static int vertex_is_integral(__isl_keep isl_basic_set *vertex) -{ - unsigned nvar; - unsigned nparam; - int i; - - nvar = isl_basic_set_dim(vertex, isl_dim_set); - nparam = isl_basic_set_dim(vertex, isl_dim_param); - for (i = 0; i < nvar; ++i) { - int r = nvar - 1 - i; - if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) && - !isl_int_is_negone(vertex->eq[r][1 + nparam + i])) - return 0; - } - - return 1; -} - -static __isl_give isl_qpolynomial *vertex_coordinate( - __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *dim) -{ - unsigned nvar; - unsigned nparam; - int r; - isl_int denom; - isl_qpolynomial *v; - - nvar = isl_basic_set_dim(vertex, isl_dim_set); - nparam = isl_basic_set_dim(vertex, isl_dim_param); - r = nvar - 1 - i; - - isl_int_init(denom); - isl_int_set(denom, vertex->eq[r][1 + nparam + i]); - isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error); - - if (isl_int_is_pos(denom)) - isl_seq_neg(vertex->eq[r], vertex->eq[r], - 1 + isl_basic_set_total_dim(vertex)); - else - isl_int_neg(denom, denom); - - v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom); - isl_int_clear(denom); - - return v; -error: - isl_space_free(dim); - isl_int_clear(denom); - return NULL; -} - -/* Check whether the bound associated to the selection "k" is tight, - * which is the case if we select exactly one vertex and if that vertex - * is integral for all values of the parameters. - */ -static int is_tight(int *k, int n, int d, isl_cell *cell) -{ - int i; - - for (i = 0; i < n; ++i) { - int v; - if (k[i] != d) { - if (k[i]) - return 0; - continue; - } - v = cell->ids[n - 1 - i]; - return vertex_is_integral(cell->vertices->v[v].vertex); - } - - return 0; -} - -static void add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom, - int *k, int n, int d, struct bernstein_data *data) -{ - isl_qpolynomial_fold *fold; - - fold = isl_qpolynomial_fold_alloc(data->type, b); - - if (data->check_tight && is_tight(k, n, d, data->cell)) - data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom, - data->fold_tight, fold); - else - data->fold = isl_qpolynomial_fold_fold_on_domain(dom, - data->fold, fold); -} - -/* Extract the coefficients of the Bernstein base polynomials and store - * them in data->fold and data->fold_tight. - * - * In particular, the coefficient of each monomial - * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding - * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]! - * - * c[i] contains the coefficient of the selected powers of the first i+1 vars. - * multinom[i] contains the partial multinomial coefficient. - */ -static void extract_coefficients(isl_qpolynomial *poly, - __isl_keep isl_set *dom, struct bernstein_data *data) -{ - int i; - int d; - int n; - isl_ctx *ctx; - isl_qpolynomial **c = NULL; - int *k = NULL; - int *left = NULL; - isl_vec *multinom = NULL; - - if (!poly) - return; - - ctx = isl_qpolynomial_get_ctx(poly); - n = isl_qpolynomial_dim(poly, isl_dim_in); - d = isl_qpolynomial_degree(poly); - isl_assert(ctx, n >= 2, return); - - c = isl_calloc_array(ctx, isl_qpolynomial *, n); - k = isl_alloc_array(ctx, int, n); - left = isl_alloc_array(ctx, int, n); - multinom = isl_vec_alloc(ctx, n); - if (!c || !k || !left || !multinom) - goto error; - - isl_int_set_si(multinom->el[0], 1); - for (k[0] = d; k[0] >= 0; --k[0]) { - int i = 1; - isl_qpolynomial_free(c[0]); - c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]); - left[0] = d - k[0]; - k[1] = -1; - isl_int_set(multinom->el[1], multinom->el[0]); - while (i > 0) { - if (i == n - 1) { - int j; - isl_space *dim; - isl_qpolynomial *b; - isl_qpolynomial *f; - for (j = 2; j <= left[i - 1]; ++j) - isl_int_divexact_ui(multinom->el[i], - multinom->el[i], j); - b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in, - n - 1 - i, left[i - 1]); - b = isl_qpolynomial_project_domain_on_params(b); - dim = isl_qpolynomial_get_domain_space(b); - f = isl_qpolynomial_rat_cst_on_domain(dim, ctx->one, - multinom->el[i]); - b = isl_qpolynomial_mul(b, f); - k[n - 1] = left[n - 2]; - add_fold(b, dom, k, n, d, data); - --i; - continue; - } - if (k[i] >= left[i - 1]) { - --i; - continue; - } - ++k[i]; - if (k[i]) - isl_int_divexact_ui(multinom->el[i], - multinom->el[i], k[i]); - isl_qpolynomial_free(c[i]); - c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in, - n - 1 - i, k[i]); - left[i] = left[i - 1] - k[i]; - k[i + 1] = -1; - isl_int_set(multinom->el[i + 1], multinom->el[i]); - ++i; - } - isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]); - } - - for (i = 0; i < n; ++i) - isl_qpolynomial_free(c[i]); - - isl_vec_free(multinom); - free(left); - free(k); - free(c); - return; -error: - isl_vec_free(multinom); - free(left); - free(k); - if (c) - for (i = 0; i < n; ++i) - isl_qpolynomial_free(c[i]); - free(c); - return; -} - -/* Perform bernstein expansion on the parametric vertices that are active - * on "cell". - * - * data->poly has been homogenized in the calling function. - * - * We plug in the barycentric coordinates for the set variables - * - * \vec x = \sum_i \alpha_i v_i(\vec p) - * - * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension. - * Next, we extract the coefficients of the Bernstein base polynomials. - */ -static int bernstein_coefficients_cell(__isl_take isl_cell *cell, void *user) -{ - int i, j; - struct bernstein_data *data = (struct bernstein_data *)user; - isl_space *dim_param; - isl_space *dim_dst; - isl_qpolynomial *poly = data->poly; - unsigned nvar; - int n_vertices; - isl_qpolynomial **subs; - isl_pw_qpolynomial_fold *pwf; - isl_set *dom; - isl_ctx *ctx; - - if (!poly) - goto error; - - nvar = isl_qpolynomial_dim(poly, isl_dim_in) - 1; - n_vertices = cell->n_vertices; - - ctx = isl_qpolynomial_get_ctx(poly); - if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate) - return isl_cell_foreach_simplex(cell, - &bernstein_coefficients_cell, user); - - subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar); - if (!subs) - goto error; - - dim_param = isl_basic_set_get_space(cell->dom); - dim_dst = isl_qpolynomial_get_domain_space(poly); - dim_dst = isl_space_add_dims(dim_dst, isl_dim_set, n_vertices); - - for (i = 0; i < 1 + nvar; ++i) - subs[i] = isl_qpolynomial_zero_on_domain(isl_space_copy(dim_dst)); - - for (i = 0; i < n_vertices; ++i) { - isl_qpolynomial *c; - c = isl_qpolynomial_var_on_domain(isl_space_copy(dim_dst), isl_dim_set, - 1 + nvar + i); - for (j = 0; j < nvar; ++j) { - int k = cell->ids[i]; - isl_qpolynomial *v; - v = vertex_coordinate(cell->vertices->v[k].vertex, j, - isl_space_copy(dim_param)); - v = isl_qpolynomial_add_dims(v, isl_dim_in, - 1 + nvar + n_vertices); - v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c)); - subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v); - } - subs[0] = isl_qpolynomial_add(subs[0], c); - } - isl_space_free(dim_dst); - - poly = isl_qpolynomial_copy(poly); - - poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices); - poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs); - poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar); - - data->cell = cell; - dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom)); - data->fold = isl_qpolynomial_fold_empty(data->type, isl_space_copy(dim_param)); - data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param); - extract_coefficients(poly, dom, data); - - pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom), - data->fold); - data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf); - pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight); - data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf); - - isl_qpolynomial_free(poly); - isl_cell_free(cell); - for (i = 0; i < 1 + nvar; ++i) - isl_qpolynomial_free(subs[i]); - free(subs); - return 0; -error: - isl_cell_free(cell); - return -1; -} - -/* Base case of applying bernstein expansion. - * - * We compute the chamber decomposition of the parametric polytope "bset" - * and then perform bernstein expansion on the parametric vertices - * that are active on each chamber. - */ -static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base( - __isl_take isl_basic_set *bset, - __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) -{ - unsigned nvar; - isl_space *dim; - isl_pw_qpolynomial_fold *pwf; - isl_vertices *vertices; - int covers; - - nvar = isl_basic_set_dim(bset, isl_dim_set); - if (nvar == 0) { - isl_set *dom; - isl_qpolynomial_fold *fold; - - fold = isl_qpolynomial_fold_alloc(data->type, poly); - dom = isl_set_from_basic_set(bset); - if (tight) - *tight = 1; - pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold); - return isl_pw_qpolynomial_fold_project_domain_on_params(pwf); - } - - if (isl_qpolynomial_is_zero(poly)) { - isl_set *dom; - isl_qpolynomial_fold *fold; - fold = isl_qpolynomial_fold_alloc(data->type, poly); - dom = isl_set_from_basic_set(bset); - pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold); - if (tight) - *tight = 1; - return isl_pw_qpolynomial_fold_project_domain_on_params(pwf); - } - - dim = isl_basic_set_get_space(bset); - dim = isl_space_params(dim); - dim = isl_space_from_domain(dim); - dim = isl_space_add_dims(dim, isl_dim_set, 1); - data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(dim), data->type); - data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type); - data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly)); - vertices = isl_basic_set_compute_vertices(bset); - isl_vertices_foreach_disjoint_cell(vertices, - &bernstein_coefficients_cell, data); - isl_vertices_free(vertices); - isl_qpolynomial_free(data->poly); - - isl_basic_set_free(bset); - isl_qpolynomial_free(poly); - - covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf); - if (covers < 0) - goto error; - - if (tight) - *tight = covers; - - if (covers) { - isl_pw_qpolynomial_fold_free(data->pwf); - return data->pwf_tight; - } - - data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight); - - return data->pwf; -error: - isl_pw_qpolynomial_fold_free(data->pwf_tight); - isl_pw_qpolynomial_fold_free(data->pwf); - return NULL; -} - -/* Apply bernstein expansion recursively by working in on len[i] - * set variables at a time, with i ranging from n_group - 1 to 0. - */ -static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive( - __isl_take isl_pw_qpolynomial *pwqp, - int n_group, int *len, struct bernstein_data *data, int *tight) -{ - int i; - unsigned nparam; - unsigned nvar; - isl_pw_qpolynomial_fold *pwf; - - if (!pwqp) - return NULL; - - nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param); - nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in); - - pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam, - isl_dim_in, 0, nvar - len[n_group - 1]); - pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight); - - for (i = n_group - 2; i >= 0; --i) { - nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param); - pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0, - isl_dim_param, nparam - len[i], len[i]); - if (tight && !*tight) - tight = NULL; - pwf = isl_pw_qpolynomial_fold_bound(pwf, tight); - } - - return pwf; -} - -static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors( - __isl_take isl_basic_set *bset, - __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) -{ - isl_factorizer *f; - isl_set *set; - isl_pw_qpolynomial *pwqp; - isl_pw_qpolynomial_fold *pwf; - - f = isl_basic_set_factorizer(bset); - if (!f) - goto error; - if (f->n_group == 0) { - isl_factorizer_free(f); - return bernstein_coefficients_base(bset, poly, data, tight); - } - - set = isl_set_from_basic_set(bset); - pwqp = isl_pw_qpolynomial_alloc(set, poly); - pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph)); - - pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data, - tight); - - isl_factorizer_free(f); - - return pwf; -error: - isl_basic_set_free(bset); - isl_qpolynomial_free(poly); - return NULL; -} - -static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive( - __isl_take isl_basic_set *bset, - __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight) -{ - int i; - int *len; - unsigned nvar; - isl_pw_qpolynomial_fold *pwf; - isl_set *set; - isl_pw_qpolynomial *pwqp; - - if (!bset || !poly) - goto error; - - nvar = isl_basic_set_dim(bset, isl_dim_set); - - len = isl_alloc_array(bset->ctx, int, nvar); - if (!len) - goto error; - - for (i = 0; i < nvar; ++i) - len[i] = 1; - - set = isl_set_from_basic_set(bset); - pwqp = isl_pw_qpolynomial_alloc(set, poly); - - pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight); - - free(len); - - return pwf; -error: - isl_basic_set_free(bset); - isl_qpolynomial_free(poly); - return NULL; -} - -/* Compute a bound on the polynomial defined over the parametric polytope - * using bernstein expansion and store the result - * in bound->pwf and bound->pwf_tight. - * - * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if - * the polytope can be factorized and apply bernstein expansion recursively - * on the factors. - * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply - * bernstein expansion recursively on each dimension. - * Otherwise, we apply bernstein expansion on the entire polytope. - */ -int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set *bset, - __isl_take isl_qpolynomial *poly, struct isl_bound *bound) -{ - struct bernstein_data data; - isl_pw_qpolynomial_fold *pwf; - unsigned nvar; - int tight = 0; - int *tp = bound->check_tight ? &tight : NULL; - - if (!bset || !poly) - goto error; - - data.type = bound->type; - data.check_tight = bound->check_tight; - - nvar = isl_basic_set_dim(bset, isl_dim_set); - - if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS) - pwf = bernstein_coefficients_factors(bset, poly, &data, tp); - else if (nvar > 1 && - (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS)) - pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp); - else - pwf = bernstein_coefficients_base(bset, poly, &data, tp); - - if (tight) - bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf); - else - bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf); - - return 0; -error: - isl_basic_set_free(bset); - isl_qpolynomial_free(poly); - return -1; -} |