forked from OSchip/llvm-project
559 lines
15 KiB
C
559 lines
15 KiB
C
/*
|
|
* Copyright 2006-2007 Universiteit Leiden
|
|
* Copyright 2008-2009 Katholieke Universiteit Leuven
|
|
* Copyright 2010 INRIA Saclay
|
|
*
|
|
* Use of this software is governed by the MIT 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_vec_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 isl_stat 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 isl_stat_ok;
|
|
error:
|
|
isl_cell_free(cell);
|
|
return isl_stat_error;
|
|
}
|
|
|
|
/* 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);
|
|
if (isl_vertices_foreach_disjoint_cell(vertices,
|
|
&bernstein_coefficients_cell, data) < 0)
|
|
data->pwf = isl_pw_qpolynomial_fold_free(data->pwf);
|
|
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 (nvar && !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.
|
|
*/
|
|
isl_stat 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 isl_stat_ok;
|
|
error:
|
|
isl_basic_set_free(bset);
|
|
isl_qpolynomial_free(poly);
|
|
return isl_stat_error;
|
|
}
|