add isl_pw_qpolynomial_bound_range
authorSven Verdoolaege <skimo@kotnet.org>
Sat, 8 May 2010 09:43:15 +0000 (11:43 +0200)
committerSven Verdoolaege <skimo@kotnet.org>
Fri, 14 May 2010 10:59:01 +0000 (12:59 +0200)
Makefile.am
include/isl_polynomial.h
isl_range.c [new file with mode: 0644]

index ee14250..c70f33e 100644 (file)
@@ -74,6 +74,7 @@ libisl_la_SOURCES = \
        isl_polynomial.c \
        isl_printer_private.h \
        isl_printer.c \
+       isl_range.c \
        isl_sample.h \
        isl_sample.c \
        isl_scan.c \
index b2e93d8..ae6cb01 100644 (file)
@@ -263,6 +263,9 @@ __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);
+
 #if defined(__cplusplus)
 }
 #endif
diff --git a/isl_range.c b/isl_range.c
new file mode 100644 (file)
index 0000000..d0e342e
--- /dev/null
@@ -0,0 +1,627 @@
+#include <isl_constraint.h>
+#include <isl_set.h>
+#include <isl_polynomial_private.h>
+#include <isl_morph.h>
+
+struct range_data {
+       int                     *signs;
+       int                     sign;
+       int                     test_monotonicity;
+       int                     monotonicity;
+       int                     exact;
+       isl_qpolynomial         *qp;
+       isl_qpolynomial         *poly;
+       isl_pw_qpolynomial_fold *pwf;
+       isl_pw_qpolynomial_fold *pwf_exact;
+};
+
+static int propagate_on_domain(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, struct range_data *data);
+
+/* Check whether the polynomial "poly" has sign "sign" over "bset",
+ * i.e., if sign == 1, check that the lower bound on the polynomial
+ * is non-negative and if sign == -1, check that the upper bound on
+ * the polynomial is non-positive.
+ */
+static int has_sign(__isl_keep isl_basic_set *bset,
+       __isl_keep isl_qpolynomial *poly, int sign, int *signs)
+{
+       struct range_data data_m;
+       unsigned nvar;
+       unsigned nparam;
+       isl_dim *dim;
+       isl_qpolynomial *opt;
+       int r;
+
+       nparam = isl_basic_set_dim(bset, isl_dim_param);
+       nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+       bset = isl_basic_set_copy(bset);
+       poly = isl_qpolynomial_copy(poly);
+
+       bset = isl_basic_set_move_dims(bset, isl_dim_set, 0,
+                                       isl_dim_param, 0, nparam);
+       poly = isl_qpolynomial_move_dims(poly, isl_dim_set, 0,
+                                       isl_dim_param, 0, nparam);
+
+       dim = isl_qpolynomial_get_dim(poly);
+       dim = isl_dim_drop(dim, isl_dim_set, 0, isl_dim_size(dim, isl_dim_set));
+
+       data_m.test_monotonicity = 0;
+       data_m.signs = signs;
+       data_m.pwf = isl_pw_qpolynomial_fold_zero(dim);
+       data_m.sign = -sign;
+       data_m.exact = 0;
+       data_m.pwf_exact = NULL;
+
+       if (propagate_on_domain(bset, poly, &data_m) < 0)
+               goto error;
+
+       if (sign > 0)
+               opt = isl_pw_qpolynomial_fold_min(data_m.pwf);
+       else
+               opt = isl_pw_qpolynomial_fold_max(data_m.pwf);
+
+       if (!opt)
+               r = -1;
+       else if (isl_qpolynomial_is_nan(opt) ||
+                isl_qpolynomial_is_infty(opt) ||
+                isl_qpolynomial_is_neginfty(opt))
+               r = 0;
+       else
+               r = sign * isl_qpolynomial_sgn(opt) >= 0;
+
+       isl_qpolynomial_free(opt);
+
+       return r;
+error:
+       isl_pw_qpolynomial_fold_free(data_m.pwf);
+       return -1;
+}
+
+/* Return  1 if poly is monotonically increasing in the last set variable,
+ *        -1 if poly is monotonically decreasing in the last set variable,
+ *        0 if no conclusion,
+ *       -2 on error.
+ *
+ * We simply check the sign of p(x+1)-p(x)
+ */
+static int monotonicity(__isl_keep isl_basic_set *bset,
+       __isl_keep isl_qpolynomial *poly, struct range_data *data)
+{
+       isl_ctx *ctx;
+       isl_dim *dim;
+       isl_qpolynomial *sub = NULL;
+       isl_qpolynomial *diff = NULL;
+       int result = 0;
+       int s;
+       unsigned nvar;
+
+       ctx = isl_qpolynomial_get_ctx(poly);
+       dim = isl_qpolynomial_get_dim(poly);
+
+       nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+       sub = isl_qpolynomial_var(isl_dim_copy(dim), isl_dim_set, nvar - 1);
+       sub = isl_qpolynomial_add(sub,
+               isl_qpolynomial_rat_cst(dim, ctx->one, ctx->one));
+
+       diff = isl_qpolynomial_substitute(isl_qpolynomial_copy(poly),
+                       isl_dim_set, nvar - 1, 1, &sub);
+       diff = isl_qpolynomial_sub(diff, isl_qpolynomial_copy(poly));
+
+       s = has_sign(bset, diff, 1, data->signs);
+       if (s < 0)
+               goto error;
+       if (s)
+               result = 1;
+       else {
+               s = has_sign(bset, diff, -1, data->signs);
+               if (s < 0)
+                       goto error;
+               if (s)
+                       result = -1;
+       }
+
+       isl_qpolynomial_free(diff);
+       isl_qpolynomial_free(sub);
+
+       return result;
+error:
+       isl_qpolynomial_free(diff);
+       isl_qpolynomial_free(sub);
+       return -2;
+}
+
+static __isl_give isl_qpolynomial *bound2poly(__isl_take isl_constraint *bound,
+       __isl_take isl_dim *dim, unsigned pos, int sign)
+{
+       if (!bound) {
+               if (sign > 0)
+                       return isl_qpolynomial_infty(dim);
+               else
+                       return isl_qpolynomial_neginfty(dim);
+       }
+       isl_dim_free(dim);
+       return isl_qpolynomial_from_constraint(bound, isl_dim_set, pos);
+}
+
+static int bound_is_integer(__isl_take isl_constraint *bound, unsigned pos)
+{
+       isl_int c;
+       int is_int;
+
+       if (!bound)
+               return 1;
+
+       isl_int_init(c);
+       isl_constraint_get_coefficient(bound, isl_dim_set, pos, &c);
+       is_int = isl_int_is_one(c) || isl_int_is_negone(c);
+       isl_int_clear(c);
+
+       return is_int;
+}
+
+struct isl_fixed_sign_data {
+       int             *signs;
+       int             sign;
+       isl_qpolynomial *poly;
+};
+
+/* Add term "term" to data->poly if it has sign data->sign.
+ * The sign is determined based on the signs of the parameters
+ * and variables in data->signs.
+ */
+static int collect_fixed_sign_terms(__isl_take isl_term *term, void *user)
+{
+       struct isl_fixed_sign_data *data = (struct isl_fixed_sign_data *)user;
+       isl_int n, d;
+       int i;
+       int sign;
+       unsigned nparam;
+       unsigned nvar;
+
+       if (!term)
+               return -1;
+
+       nparam = isl_term_dim(term, isl_dim_param);
+       nvar = isl_term_dim(term, isl_dim_set);
+
+       isl_assert(isl_term_get_ctx(term), isl_term_dim(term, isl_dim_div) == 0,
+                       return -1);
+
+       isl_int_init(n);
+       isl_int_init(d);
+
+       isl_term_get_num(term, &n);
+       isl_term_get_den(term, &d);
+
+       sign = isl_int_sgn(n);
+       for (i = 0; i < nparam; ++i) {
+               if (data->signs[i] > 0)
+                       continue;
+               if (isl_term_get_exp(term, isl_dim_param, i) % 2)
+                       sign = -sign;
+       }
+       for (i = 0; i < nvar; ++i) {
+               if (data->signs[nparam + i] > 0)
+                       continue;
+               if (isl_term_get_exp(term, isl_dim_set, i) % 2)
+                       sign = -sign;
+       }
+
+       if (sign == data->sign) {
+               isl_qpolynomial *t = isl_qpolynomial_from_term(term);
+
+               data->poly = isl_qpolynomial_add(data->poly, t);
+       } else
+               isl_term_free(term);
+
+       isl_int_clear(n);
+       isl_int_clear(d);
+
+       return 0;
+}
+
+/* Construct and return a polynomial that consists of the terms
+ * in "poly" that have sign "sign".
+ */
+static __isl_give isl_qpolynomial *fixed_sign_terms(
+       __isl_keep isl_qpolynomial *poly, int *signs, int sign)
+{
+       struct isl_fixed_sign_data data = { signs, sign };
+       data.poly = isl_qpolynomial_zero(isl_qpolynomial_get_dim(poly));
+
+       if (isl_qpolynomial_foreach_term(poly, collect_fixed_sign_terms, &data) < 0)
+               goto error;
+
+       return data.poly;
+error:
+       isl_qpolynomial_free(data.poly);
+       return NULL;
+}
+
+/* Helper function to add a guarder polynomial to either pwf_exact or pwf,
+ * depending on whether the result has been determined to be exact.
+ */
+static int add_guarded_poly(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, struct range_data *data)
+{
+       enum isl_fold type = data->sign < 0 ? isl_fold_min : isl_fold_max;
+       isl_set *set;
+       isl_qpolynomial_fold *fold;
+       isl_pw_qpolynomial_fold *pwf;
+
+       fold = isl_qpolynomial_fold_alloc(type, poly);
+       set = isl_set_from_basic_set(bset);
+       pwf = isl_pw_qpolynomial_fold_alloc(set, fold);
+       if (data->exact)
+               data->pwf_exact = isl_pw_qpolynomial_fold_add(
+                                               data->pwf_exact, pwf);
+       else
+               data->pwf = isl_pw_qpolynomial_fold_add(data->pwf, pwf);
+
+       return 0;
+}
+
+/* Given a lower and upper bound on the final variable and constraints
+ * on the remaining variables where these bounds are active,
+ * eliminate the variable from data->poly based on these bounds.
+ * If the polynomial has been determined to be monotonic
+ * in the variable, then simply plug in the appropriate bound.
+ * If the current polynomial is exact and if this bound is integer,
+ * then the result is still exact.  In all other cases, the results
+ * may be inexact.
+ * Otherwise, plug in the largest bound (in absolute value) in
+ * the positive terms (if an upper bound is wanted) or the negative terms
+ * (if a lower bounded is wanted) and the other bound in the other terms.
+ *
+ * If all variables have been eliminated, then record the result.
+ * Ohterwise, recurse on the next variable.
+ */
+static int propagate_on_bound_pair(__isl_take isl_constraint *lower,
+       __isl_take isl_constraint *upper, __isl_take isl_basic_set *bset,
+       void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+       int save_exact = data->exact;
+       isl_qpolynomial *poly;
+       int r;
+       unsigned nvar;
+
+       nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+       if (data->monotonicity) {
+               isl_qpolynomial *sub;
+               isl_dim *dim = isl_qpolynomial_get_dim(data->poly);
+               if (data->monotonicity * data->sign > 0) {
+                       if (data->exact)
+                               data->exact = bound_is_integer(upper, nvar);
+                       sub = bound2poly(upper, dim, nvar, 1);
+                       isl_constraint_free(lower);
+               } else {
+                       if (data->exact)
+                               data->exact = bound_is_integer(lower, nvar);
+                       sub = bound2poly(lower, dim, nvar, -1);
+                       isl_constraint_free(upper);
+               }
+               poly = isl_qpolynomial_copy(data->poly);
+               poly = isl_qpolynomial_substitute(poly, isl_dim_set, nvar, 1, &sub);
+               poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, nvar, 1);
+
+               isl_qpolynomial_free(sub);
+       } else {
+               isl_qpolynomial *l, *u;
+               isl_qpolynomial *pos, *neg;
+               isl_dim *dim = isl_qpolynomial_get_dim(data->poly);
+               unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
+               int sign = data->sign * data->signs[nparam + nvar];
+
+               data->exact = 0;
+
+               u = bound2poly(upper, isl_dim_copy(dim), nvar, 1);
+               l = bound2poly(lower, dim, nvar, -1);
+
+               pos = fixed_sign_terms(data->poly, data->signs, sign);
+               neg = fixed_sign_terms(data->poly, data->signs, -sign);
+
+               pos = isl_qpolynomial_substitute(pos, isl_dim_set, nvar, 1, &u);
+               neg = isl_qpolynomial_substitute(neg, isl_dim_set, nvar, 1, &l);
+
+               poly = isl_qpolynomial_add(pos, neg);
+               poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, nvar, 1);
+
+               isl_qpolynomial_free(u);
+               isl_qpolynomial_free(l);
+       }
+
+       if (isl_basic_set_dim(bset, isl_dim_set) == 0)
+               r = add_guarded_poly(bset, poly, data);
+       else
+               r = propagate_on_domain(bset, poly, data);
+
+       data->exact = save_exact;
+
+       return r;
+}
+
+/* Recursively perform range propagation on the polynomial "poly"
+ * defined over the basic set "bset" and collect the results in "data".
+ */
+static int propagate_on_domain(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, struct range_data *data)
+{
+       isl_qpolynomial *save_poly = data->poly;
+       int save_monotonicity = data->monotonicity;
+       unsigned d;
+
+       if (!bset || !poly)
+               goto error;
+
+       d = isl_basic_set_dim(bset, isl_dim_set);
+       isl_assert(bset->ctx, d >= 1, goto error);
+
+       if (isl_qpolynomial_is_cst(poly, NULL, NULL)) {
+               bset = isl_basic_set_project_out(bset, isl_dim_set, 0, d);
+               poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, 0, d);
+               return add_guarded_poly(bset, poly, data);
+       }
+
+       if (data->test_monotonicity)
+               data->monotonicity = monotonicity(bset, poly, data);
+       else
+               data->monotonicity = 0;
+       if (data->monotonicity < -1)
+               goto error;
+
+       data->poly = poly;
+       if (isl_basic_set_foreach_bound_pair(bset, isl_dim_set, d - 1,
+                                           &propagate_on_bound_pair, data) < 0)
+               goto error;
+
+       isl_basic_set_free(bset);
+       isl_qpolynomial_free(poly);
+       data->monotonicity = save_monotonicity;
+       data->poly = save_poly;
+
+       return 0;
+error:
+       isl_basic_set_free(bset);
+       isl_qpolynomial_free(poly);
+       data->monotonicity = save_monotonicity;
+       data->poly = save_poly;
+       return -1;
+}
+
+static int basic_guarded_poly_bound(__isl_take isl_basic_set *bset, void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+       unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
+       unsigned dim = isl_basic_set_dim(bset, isl_dim_set);
+       int r;
+
+       data->signs = NULL;
+
+       data->signs = isl_alloc_array(bset->ctx, int,
+                                       isl_basic_set_dim(bset, isl_dim_all));
+
+       if (isl_basic_set_dims_get_sign(bset, isl_dim_set, 0, dim,
+                                       data->signs + nparam) < 0)
+               goto error;
+       if (isl_basic_set_dims_get_sign(bset, isl_dim_param, 0, nparam,
+                                       data->signs) < 0)
+               goto error;
+
+       r = propagate_on_domain(bset, isl_qpolynomial_copy(data->poly), data);
+
+       free(data->signs);
+
+       return r;
+error:
+       free(data->signs);
+       isl_basic_set_free(bset);
+       return -1;
+}
+
+static int compressed_guarded_poly_bound(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+       unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
+       unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
+       isl_set *set;
+
+       if (!bset)
+               goto error;
+
+       if (nvar == 0)
+               return add_guarded_poly(bset, poly, data);
+
+       set = isl_set_from_basic_set(bset);
+       set = isl_set_split_dims(set, isl_dim_param, 0, nparam);
+       set = isl_set_split_dims(set, isl_dim_set, 0, nvar);
+
+       data->poly = poly;
+
+       if (isl_set_foreach_basic_set(set, &basic_guarded_poly_bound, data) < 0)
+               goto error;
+
+       isl_set_free(set);
+       isl_qpolynomial_free(poly);
+
+       return 0;
+error:
+       isl_set_free(set);
+       isl_qpolynomial_free(poly);
+       return -1;
+}
+
+static int guarded_poly_bound(__isl_take isl_basic_set *bset,
+       __isl_take isl_qpolynomial *poly, void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+       isl_pw_qpolynomial_fold *top_pwf;
+       isl_pw_qpolynomial_fold *top_pwf_exact;
+       isl_dim *dim;
+       isl_morph *morph, *morph2;
+       unsigned orig_nvar;
+       int r;
+
+       bset = isl_basic_set_detect_equalities(bset);
+
+       if (!bset)
+               goto error;
+
+       if (bset->n_eq == 0)
+               return compressed_guarded_poly_bound(bset, poly, user);
+
+       orig_nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+       morph = isl_basic_set_variable_compression(bset, isl_dim_param);
+       bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
+
+       morph2 = isl_basic_set_parameter_compression(bset);
+       bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
+
+       morph = isl_morph_compose(morph2, morph);
+
+       morph2 = isl_basic_set_variable_compression(bset, isl_dim_set);
+       bset = isl_morph_basic_set(isl_morph_copy(morph2), bset);
+
+       morph2 = isl_morph_compose(morph2, isl_morph_copy(morph));
+       poly = isl_qpolynomial_morph(poly, morph2);
+
+       dim = isl_morph_get_ran_dim(morph);
+       dim = isl_dim_drop(dim, isl_dim_set, 0, isl_dim_size(dim, isl_dim_set));
+
+       top_pwf = data->pwf;
+       top_pwf_exact = data->pwf_exact;
+
+       data->pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
+       data->pwf_exact = isl_pw_qpolynomial_fold_zero(dim);
+
+       r = compressed_guarded_poly_bound(bset, poly, user);
+
+       morph = isl_morph_drop_dims(morph, isl_dim_set, 0, orig_nvar);
+       morph = isl_morph_inverse(morph);
+
+       data->pwf = isl_pw_qpolynomial_fold_morph(data->pwf,
+                                                       isl_morph_copy(morph));
+       data->pwf_exact = isl_pw_qpolynomial_fold_morph(data->pwf_exact, morph);
+
+       data->pwf = isl_pw_qpolynomial_fold_add(top_pwf, data->pwf);
+       data->pwf_exact = isl_pw_qpolynomial_fold_add(top_pwf_exact,
+                                                       data->pwf_exact);
+
+       return r;
+error:
+       isl_basic_set_free(bset);
+       isl_qpolynomial_free(poly);
+       return -1;
+}
+
+static int basic_guarded_bound(__isl_take isl_basic_set *bset, void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+       int r;
+
+       r = isl_qpolynomial_as_polynomial_on_domain(data->qp, bset,
+                                                   &guarded_poly_bound, user);
+       isl_basic_set_free(bset);
+       return r;
+}
+
+static int guarded_bound(__isl_take isl_set *set,
+       __isl_take isl_qpolynomial *qp, void *user)
+{
+       struct range_data *data = (struct range_data *)user;
+
+       if (!set || !qp)
+               goto error;
+
+       set = isl_set_make_disjoint(set);
+
+       data->qp = qp;
+
+       if (isl_set_foreach_basic_set(set, &basic_guarded_bound, data) < 0)
+               goto error;
+
+       isl_set_free(set);
+       isl_qpolynomial_free(qp);
+
+       return 0;
+error:
+       isl_set_free(set);
+       isl_qpolynomial_free(qp);
+       return -1;
+}
+
+/* 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 *exact)
+{
+       isl_dim *dim;
+       isl_pw_qpolynomial_fold *pwf;
+       unsigned nvar;
+       unsigned nparam;
+       struct range_data data;
+       int covers;
+
+       if (!pwqp)
+               return NULL;
+
+       dim = isl_pw_qpolynomial_get_dim(pwqp);
+       nvar = isl_dim_size(dim, isl_dim_set);
+
+       if (isl_pw_qpolynomial_is_zero(pwqp)) {
+               isl_pw_qpolynomial_free(pwqp);
+               dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
+               return isl_pw_qpolynomial_fold_zero(dim);
+       }
+
+       if (nvar == 0) {
+               isl_dim_free(dim);
+               return isl_pw_qpolynomial_fold_from_pw_qpolynomial(type, pwqp);
+       }
+
+       dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
+
+       nparam = isl_dim_size(dim, isl_dim_param);
+       data.pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
+       data.pwf_exact = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
+       if (type == isl_fold_min)
+               data.sign = -1;
+       else
+               data.sign = 1;
+       data.test_monotonicity = 1;
+       data.exact = !!exact;
+
+       if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp, guarded_bound, &data))
+               goto error;
+
+       covers = isl_pw_qpolynomial_fold_covers(data.pwf_exact, data.pwf);
+       if (covers < 0)
+               goto error;
+
+       if (exact)
+               *exact = covers;
+
+       isl_dim_free(dim);
+       isl_pw_qpolynomial_free(pwqp);
+
+       if (covers) {
+               isl_pw_qpolynomial_fold_free(data.pwf);
+               return data.pwf_exact;
+       }
+
+       data.pwf = isl_pw_qpolynomial_fold_add(data.pwf, data.pwf_exact);
+
+       return data.pwf;
+error:
+       isl_pw_qpolynomial_fold_free(data.pwf);
+       isl_dim_free(dim);
+       isl_pw_qpolynomial_free(pwqp);
+       return NULL;
+}