isl_convex_hull.c: explicitly mark wrapping lp problem as rational

[isl.git] / isl_range.c

#include <isl/constraint.h>
#include <isl/set.h>
#include <isl_polynomial_private.h>
#include <isl_morph.h>
#include <isl_range.h>

struct range_data {
        struct isl_bound        *bound;
        int                     *signs;
        int                     sign;
        int                     test_monotonicity;
        int                     monotonicity;
        int                     tight;
        isl_qpolynomial         *poly;
        isl_pw_qpolynomial_fold *pwf;
        isl_pw_qpolynomial_fold *pwf_tight;
};

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;
        enum isl_fold type;

        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.sign = -sign;
        type = data_m.sign < 0 ? isl_fold_min : isl_fold_max;
        data_m.pwf = isl_pw_qpolynomial_fold_zero(dim, type);
        data_m.tight = 0;
        data_m.pwf_tight = 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.  The integer divisions, if
 * any, are assumed to be non-negative.
 */
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;
        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_int_init(n);

        isl_term_get_num(term, &n);

        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);

        return 0;
}

/* Construct and return a polynomial that consists of the terms
 * in "poly" that have sign "sign".  The integer divisions, if
 * any, are assumed to be non-negative.
 */
__isl_give isl_qpolynomial *isl_qpolynomial_terms_of_sign(
        __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 guarded polynomial to either pwf_tight or pwf,
 * depending on whether the result has been determined to be tight.
 */
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(type, set, fold);
        if (data->tight)
                data->pwf_tight = isl_pw_qpolynomial_fold_fold(
                                                data->pwf_tight, pwf);
        else
                data->pwf = isl_pw_qpolynomial_fold_fold(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 tight and if this bound is integer,
 * then the result is still tight.  In all other cases, the results
 * may not be tight.
 * 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_tight = data->tight;
        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->tight)
                                data->tight = bound_is_integer(upper, nvar);
                        sub = bound2poly(upper, dim, nvar, 1);
                        isl_constraint_free(lower);
                } else {
                        if (data->tight)
                                data->tight = 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->tight = 0;

                u = bound2poly(upper, isl_dim_copy(dim), nvar, 1);
                l = bound2poly(lower, dim, nvar, -1);

                pos = isl_qpolynomial_terms_of_sign(data->poly, data->signs, sign);
                neg = isl_qpolynomial_terms_of_sign(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->tight = save_tight;

        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_ctx *ctx;
        isl_qpolynomial *save_poly = data->poly;
        int save_monotonicity = data->monotonicity;
        unsigned d;

        if (!bset || !poly)
                goto error;

        ctx = isl_basic_set_get_ctx(bset);
        d = isl_basic_set_dim(bset, isl_dim_set);
        isl_assert(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;
        isl_ctx *ctx;
        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;

        ctx = isl_basic_set_get_ctx(bset);
        data->signs = isl_alloc_array(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 qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
        __isl_take isl_qpolynomial *poly, struct range_data *data)
{
        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;

        data->test_monotonicity = 1;
        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;
}

int isl_qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
        __isl_take isl_qpolynomial *poly, struct isl_bound *bound)
{
        struct range_data data;
        int r;

        data.pwf = bound->pwf;
        data.pwf_tight = bound->pwf_tight;
        data.tight = bound->check_tight;
        if (bound->type == isl_fold_min)
                data.sign = -1;
        else
                data.sign = 1;

        r = qpolynomial_bound_on_domain_range(bset, poly, &data);

        bound->pwf = data.pwf;
        bound->pwf_tight = data.pwf_tight;

        return r;
}