optionally (and by default) use bernstein expansion to compute bounds
authorSven Verdoolaege <skimo@kotnet.org>
Mon, 7 Jun 2010 12:43:26 +0000 (14:43 +0200)
committerSven Verdoolaege <skimo@kotnet.org>
Sat, 12 Jun 2010 15:34:01 +0000 (17:34 +0200)
12 files changed:
Makefile.am
bound.c
bound_test.sh
doc/user.pod
include/isl_options.h
include/isl_polynomial.h
isl_bernstein.c [new file with mode: 0644]
isl_bernstein.h [new file with mode: 0644]
isl_bound.c
isl_bound.h
isl_options.c
isl_range.c

index 2a4ab92..05887b1 100644 (file)
@@ -35,6 +35,8 @@ libisl_la_SOURCES = \
        isl_arg.c \
        isl_basis_reduction.h \
        basis_reduction_tab.c \
+       isl_bernstein.c \
+       isl_bernstein.h \
        isl_blk.c \
        isl_bound.c \
        isl_bound.h \
diff --git a/bound.c b/bound.c
index f3ba127..add58dc 100644 (file)
--- a/bound.c
+++ b/bound.c
@@ -243,7 +243,7 @@ int main(int argc, char **argv)
        if (options->verify)
                copy = isl_pw_qpolynomial_copy(pwqp);
 
-       pwf = isl_pw_qpolynomial_bound_range(pwqp, isl_fold_max, &exact);
+       pwf = isl_pw_qpolynomial_bound(pwqp, isl_fold_max, &exact);
        pwf = isl_pw_qpolynomial_fold_coalesce(pwf);
 
        if (options->verify) {
index 28f27e7..f26dd1b 100755 (executable)
@@ -28,5 +28,6 @@ BOUND_TESTS="\
 
 for i in $BOUND_TESTS; do
        echo $i;
-       ./isl_bound$EXEEXT -T < $srcdir/test_inputs/$i || exit
+       ./isl_bound$EXEEXT -T --bound=bernstein < $srcdir/test_inputs/$i || exit
+       ./isl_bound$EXEEXT -T --bound=range < $srcdir/test_inputs/$i || exit
 done
index 5fe3d3d..963a6ee 100644 (file)
@@ -1452,6 +1452,97 @@ the cells in the domain of the input piecewise quasipolynomial.
 In future, the operation will also exploit the context
 to simplify the quasipolynomials associated to each cell.
 
+=head2 Bounds on Piecewise Quasipolynomials and Piecewise Quasipolynomial Reductions
+
+A piecewise quasipolynomial reduction is a piecewise
+reduction (or fold) of quasipolynomials.
+In particular, the reduction can be maximum or a minimum.
+The objects are mainly used to represent the result of
+an upper or lower bound on a quasipolynomial over its domain,
+i.e., as the result of the following function.
+
+       __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound(
+               __isl_take isl_pw_qpolynomial *pwqp,
+               enum isl_fold type, int *tight);
+
+The C<type> argument may be either C<isl_fold_min> or C<isl_fold_max>.
+If C<tight> is not C<NULL>, then C<*tight> is set to C<1>
+is the returned bound is known be tight, i.e., for each value
+of the parameters there is at least
+one element in the domain that reaches the bound.
+
+A (piecewise) quasipolynomial reduction can be copied or freed using the
+following functions.
+
+       __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
+               __isl_keep isl_qpolynomial_fold *fold);
+       __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_copy(
+               __isl_keep isl_pw_qpolynomial_fold *pwf);
+       void isl_qpolynomial_fold_free(
+               __isl_take isl_qpolynomial_fold *fold);
+       void isl_pw_qpolynomial_fold_free(
+               __isl_take isl_pw_qpolynomial_fold *pwf);
+
+=head3 Printing Piecewise Quasipolynomial Reductions
+
+Piecewise quasipolynomial reductions can be printed
+using the following function.
+
+       __isl_give isl_printer *isl_printer_print_pw_qpolynomial_fold(
+               __isl_take isl_printer *p,
+               __isl_keep isl_pw_qpolynomial_fold *pwf);
+
+The output format of the printer
+needs to be set to either C<ISL_FORMAT_ISL> or C<ISL_FORMAT_C>.
+
+=head3 Inspecting (Piecewise) Quasipolynomial Reductions
+
+To iterate over the cells in a piecewise quasipolynomial reduction,
+use either of the following two functions
+
+       int isl_pw_qpolynomial_fold_foreach_piece(
+               __isl_keep isl_pw_qpolynomial_fold *pwf,
+               int (*fn)(__isl_take isl_set *set,
+                         __isl_take isl_qpolynomial_fold *fold,
+                         void *user), void *user);
+       int isl_pw_qpolynomial_fold_foreach_lifted_piece(
+               __isl_keep isl_pw_qpolynomial_fold *pwf,
+               int (*fn)(__isl_take isl_set *set,
+                         __isl_take isl_qpolynomial_fold *fold,
+                         void *user), void *user);
+
+See L<Inspecting (Piecewise) Quasipolynomials> for an explanation
+of the difference between these two functions.
+
+To iterate over all quasipolynomials in a reduction, use
+
+       int isl_qpolynomial_fold_foreach_qpolynomial(
+               __isl_keep isl_qpolynomial_fold *fold,
+               int (*fn)(__isl_take isl_qpolynomial *qp,
+                         void *user), void *user);
+
+=head3 Operations on Piecewise Quasipolynomial Reductions
+
+       __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_add(
+               __isl_take isl_pw_qpolynomial_fold *pwf1,
+               __isl_take isl_pw_qpolynomial_fold *pwf2);
+
+       __isl_give isl_qpolynomial *isl_pw_qpolynomial_fold_eval(
+               __isl_take isl_pw_qpolynomial_fold *pwf,
+               __isl_take isl_point *pnt);
+
+       __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_coalesce(
+               __isl_take isl_pw_qpolynomial_fold *pwf);
+
+       __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_gist(
+               __isl_take isl_pw_qpolynomial_fold *pwf,
+               __isl_take isl_set *context);
+
+The gist operation applies the gist operation to each of
+the cells in the domain of the input piecewise quasipolynomial reduction.
+In future, the operation will also exploit the context
+to simplify the quasipolynomial reductions associated to each cell.
+
 =head2 Dependence Analysis
 
 C<isl> contains specialized functionality for performing
index b918d03..31494ba 100644 (file)
@@ -42,6 +42,14 @@ struct isl_options {
        #define                 ISL_CLOSURE_ISL         0
        #define                 ISL_CLOSURE_OMEGA       1
        unsigned                closure;
+
+       #define                 ISL_BOUND_BERNSTEIN     0
+       #define                 ISL_BOUND_RANGE         1
+       int                     bound;
+
+       #define                 ISL_BERNSTEIN_FACTORS   1
+       #define                 ISL_BERNSTEIN_INTERVALS 2
+       int                     bernstein_recurse;
 };
 
 ISL_ARG_DECL(isl_options, struct isl_options, isl_options_arg)
index e612750..4ad8486 100644 (file)
@@ -303,8 +303,10 @@ __isl_give isl_qpolynomial *isl_pw_qpolynomial_fold_max(
 __isl_give isl_qpolynomial *isl_pw_qpolynomial_fold_min(
        __isl_take isl_pw_qpolynomial_fold *pwf);
 
-__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound_range(
-       __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type, int *exact);
+__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound(
+       __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type, int *tight);
+__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_bound(
+       __isl_take isl_pw_qpolynomial_fold *pwf, int *tight);
 
 #if defined(__cplusplus)
 }
diff --git a/isl_bernstein.c b/isl_bernstein.c
new file mode 100644 (file)
index 0000000..868bb3c
--- /dev/null
@@ -0,0 +1,540 @@
+/*
+ * 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_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_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_dim *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_dim_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, j;
+
+       for (i = 0; i < n; ++i) {
+               int v;
+               if (k[i] != d) {
+                       if (k[i])
+                               return 0;
+                       continue;
+               }
+               v = cell->vertices->c[cell->id].vertices[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_set);
+       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_set, 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_dim *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_set,
+                                       n - 1 - i, left[i - 1]);
+                               b = isl_qpolynomial_drop_dims(b, isl_dim_set,
+                                                               0, n);
+                               dim = isl_qpolynomial_get_dim(b);
+                               f = isl_qpolynomial_rat_cst(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_set,
+                                       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_dim *dim_param;
+       isl_dim *dim_dst;
+       isl_qpolynomial *poly = data->poly;
+       unsigned nvar;
+       int n_vertices;
+       isl_qpolynomial **subs;
+       isl_pw_qpolynomial_fold *pwf;
+       isl_set *dom;
+
+       nvar = isl_qpolynomial_dim(poly, isl_dim_set) - 1;
+       n_vertices = cell->vertices->c[cell->id].n_vertices;
+
+       subs = isl_alloc_array(data->poly->dim->ctx, isl_qpolynomial *,
+                               1 + nvar);
+       if (!subs)
+               goto error;
+
+       dim_param = isl_basic_set_get_dim(cell->dom);
+       dim_dst = isl_qpolynomial_get_dim(poly);
+       dim_dst = isl_dim_add(dim_dst, isl_dim_set, n_vertices);
+
+       for (i = 0; i < 1 + nvar; ++i)
+               subs[i] = isl_qpolynomial_zero(isl_dim_copy(dim_dst));
+
+       for (i = 0; i < n_vertices; ++i) {
+               isl_qpolynomial *c;
+               c = isl_qpolynomial_var(isl_dim_copy(dim_dst), isl_dim_set,
+                                       1 + nvar + i);
+               for (j = 0; j < nvar; ++j) {
+                       int k = cell->vertices->c[cell->id].vertices[i];
+                       isl_qpolynomial *v;
+                       v = vertex_coordinate(cell->vertices->v[k].vertex, j,
+                                               isl_dim_copy(dim_param));
+                       v = isl_qpolynomial_add_dims(v, isl_dim_set,
+                                                       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_dim_free(dim_dst);
+
+       poly = isl_qpolynomial_copy(poly);
+
+       poly = isl_qpolynomial_add_dims(poly, isl_dim_set, n_vertices);
+       poly = isl_qpolynomial_substitute(poly, isl_dim_set, 0, 1 + nvar, subs);
+       poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, 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_dim_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(isl_set_copy(dom), data->fold);
+       data->pwf = isl_pw_qpolynomial_fold_add(data->pwf, pwf);
+       pwf = isl_pw_qpolynomial_fold_alloc(dom, data->fold_tight);
+       data->pwf_tight = isl_pw_qpolynomial_fold_add(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_dim *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;
+               return isl_pw_qpolynomial_fold_alloc(dom, fold);
+       }
+
+       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(dom, fold);
+               if (tight)
+                       *tight = 1;
+               return isl_pw_qpolynomial_fold_drop_dims(pwf,
+                                                           isl_dim_set, 0, nvar);
+       }
+
+       dim = isl_basic_set_get_dim(bset);
+       dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
+       data->pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
+       data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim);
+       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_add(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_set);
+
+       pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,
+                                       isl_dim_set, 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_set, 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(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_add(bound->pwf_tight, pwf);
+       else
+               bound->pwf = isl_pw_qpolynomial_fold_add(bound->pwf, pwf);
+
+       return 0;
+error:
+       isl_basic_set_free(bset);
+       isl_qpolynomial_free(poly);
+       return -1;
+}
diff --git a/isl_bernstein.h b/isl_bernstein.h
new file mode 100644 (file)
index 0000000..7694b04
--- /dev/null
@@ -0,0 +1,4 @@
+#include <isl_bound.h>
+
+int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, struct isl_bound *bound);
index 40e4d5e..7660d04 100644 (file)
@@ -9,14 +9,40 @@
  */
 
 #include <isl_bound.h>
+#include <isl_bernstein.h>
 #include <isl_range.h>
 #include <isl_polynomial_private.h>
 
+/* Compute a bound on the polynomial defined over the parametric polytope
+ * using either range propagation or bernstein expansion and
+ * store the result in bound->pwf and bound->pwf_tight.
+ * Since bernstein expansion requires bounded domains, we apply
+ * range propagation on unbounded domains.  Otherwise, we respect the choice
+ * of the user.
+ */
 static int compressed_guarded_poly_bound(__isl_take isl_basic_set *bset,
        __isl_take isl_qpolynomial *poly, void *user)
 {
        struct isl_bound *bound = (struct isl_bound *)user;
-       return isl_qpolynomial_bound_on_domain_range(bset, poly, bound);
+       int bounded;
+
+       if (!bset || !poly)
+               goto error;
+
+       if (bset->ctx->opt->bound == ISL_BOUND_RANGE)
+               return isl_qpolynomial_bound_on_domain_range(bset, poly, bound);
+
+       bounded = isl_basic_set_is_bounded(bset);
+       if (bounded < 0)
+               goto error;
+       if (bounded)
+               return isl_qpolynomial_bound_on_domain_bernstein(bset, poly, bound);
+       else
+               return isl_qpolynomial_bound_on_domain_range(bset, poly, bound);
+error:
+       isl_basic_set_free(bset);
+       isl_qpolynomial_free(poly);
+       return -1;
 }
 
 static int guarded_poly_bound(__isl_take isl_basic_set *bset,
index 07a805b..b02a0b3 100644 (file)
@@ -15,7 +15,4 @@ struct isl_bound {
        isl_pw_qpolynomial_fold *pwf_tight;
 };
 
-__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound(
-       __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type, int *tight);
-
 #endif
index 2f8f0b5..92890d8 100644 (file)
@@ -57,6 +57,23 @@ struct isl_arg_choice isl_closure_choice[] = {
        {0}
 };
 
+static struct isl_arg_choice bound[] = {
+       {"bernstein",   ISL_BOUND_BERNSTEIN},
+       {"range",       ISL_BOUND_RANGE},
+       {0}
+};
+
+static struct isl_arg_flags bernstein_recurse[] = {
+       {"none",        ISL_BERNSTEIN_FACTORS | ISL_BERNSTEIN_INTERVALS, 0},
+       {"factors",     ISL_BERNSTEIN_FACTORS | ISL_BERNSTEIN_INTERVALS,
+                       ISL_BERNSTEIN_FACTORS},
+       {"intervals",   ISL_BERNSTEIN_FACTORS | ISL_BERNSTEIN_INTERVALS,
+                       ISL_BERNSTEIN_INTERVALS},
+       {"full",        ISL_BERNSTEIN_FACTORS | ISL_BERNSTEIN_INTERVALS,
+                       ISL_BERNSTEIN_FACTORS | ISL_BERNSTEIN_INTERVALS},
+       {0}
+};
+
 struct isl_arg isl_options_arg[] = {
 ISL_ARG_CHOICE(struct isl_options, lp_solver, 0, "lp-solver", \
        isl_lp_solver_choice,   ISL_LP_TAB, "lp solver to use")
@@ -75,6 +92,10 @@ ISL_ARG_CHOICE(struct isl_options, closure, 0, "closure", \
        "closure operation to use")
 ISL_ARG_BOOL(struct isl_options, gbr_only_first, 0, "gbr-only-first", 0,
        "only perform basis reduction in first direction")
+ISL_ARG_CHOICE(struct isl_options, bound, 0, "bound", bound,
+       ISL_BOUND_BERNSTEIN, "algorithm to use for computing bounds")
+ISL_ARG_FLAGS(struct isl_options, bernstein_recurse, 0,
+       "bernstein-recurse", bernstein_recurse, ISL_BERNSTEIN_FACTORS, NULL)
 ISL_ARG_END
 };
 
index 457a2df..140ca50 100644 (file)
@@ -478,12 +478,3 @@ int isl_qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
 
        return r;
 }
-
-/* Compute a lower or upper bound (depending on "type") on the given
- * piecewise step-polynomial using range propagation.
- */
-__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound_range(
-       __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type, int *tight)
-{
-       return isl_pw_qpolynomial_bound(pwqp, type, tight);
-}