--- /dev/null
+#include <assert.h>
+#include "isl_tab.h"
+
+struct tab_lp {
+ struct isl_ctx *ctx;
+ struct isl_vec *row;
+ struct isl_tab *tab;
+ struct isl_tab_undo **stack;
+ isl_int *obj;
+ isl_int opt;
+ isl_int opt_denom;
+ int neq;
+ unsigned dim;
+ int n_ineq;
+};
+
+static struct tab_lp *init_lp(struct isl_basic_set *bset);
+static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
+static int solve_lp(struct tab_lp *lp);
+static void get_obj_val(struct tab_lp* lp, mpq_t *F);
+static void delete_lp(struct tab_lp *lp);
+static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
+static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha);
+static void del_lp_row(struct tab_lp *lp);
+
+#define GBR_LP struct tab_lp
+#define GBR_type mpq_t
+#define GBR_init(v) mpq_init(v)
+#define GBR_clear(v) mpq_clear(v)
+#define GBR_set(a,b) mpq_set(a,b)
+#define GBR_set_ui(a,b) mpq_set_ui(a,b,1)
+#define GBR_mul(a,b,c) mpq_mul(a,b,c)
+#define GBR_lt(a,b) (mpq_cmp(a,b) < 0)
+#define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
+#define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
+#define GBR_lp_init(P) init_lp(P)
+#define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim)
+#define GBR_lp_solve(lp) solve_lp(lp)
+#define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F)
+#define GBR_lp_delete(lp) delete_lp(lp)
+#define GBR_lp_next_row(lp) lp->neq
+#define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim)
+#define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha)
+#define GBR_lp_del_row(lp) del_lp_row(lp);
+#include "basis_reduction_templ.c"
+
+/* Set up a tableau for the Cartesian product of bset with itself.
+ * This could be optimized by first setting up a tableau for bset
+ * and then performing the Cartesian product on the tableau.
+ */
+static struct isl_tab *gbr_tab(struct isl_basic_set *bset,
+ struct isl_vec *row)
+{
+ int i, j;
+ unsigned dim;
+ struct isl_tab *tab;
+
+ if (!bset || !row)
+ return NULL;
+
+ dim = isl_basic_set_total_dim(bset);
+ tab = isl_tab_alloc(bset->ctx, 2 * bset->n_ineq + dim + 1, 2 * dim);
+
+ for (i = 0; i < 2; ++i) {
+ isl_seq_clr(row->el + 1 + (1 - i) * dim, dim);
+ for (j = 0; j < bset->n_ineq; ++j) {
+ isl_int_set(row->el[0], bset->ineq[j][0]);
+ isl_seq_cpy(row->el + 1 + i * dim,
+ bset->ineq[j] + 1, dim);
+ tab = isl_tab_add_ineq(bset->ctx, tab, row->el);
+ if (!tab || tab->empty)
+ return tab;
+ }
+ }
+
+ return tab;
+}
+
+static struct tab_lp *init_lp(struct isl_basic_set *bset)
+{
+ struct tab_lp *lp = NULL;
+
+ if (!bset)
+ return NULL;
+
+ isl_assert(bset->ctx, bset->n_eq == 0, return NULL);
+
+ lp = isl_calloc_type(bset->ctx, struct tab_lp);
+ if (!lp)
+ return NULL;
+
+ isl_int_init(lp->opt);
+ isl_int_init(lp->opt_denom);
+
+ lp->dim = isl_basic_set_total_dim(bset);
+ lp->n_ineq = bset->n_ineq;
+
+ lp->ctx = bset->ctx;
+ isl_ctx_ref(lp->ctx);
+
+ lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
+
+ lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);
+ if (!lp->row)
+ goto error;
+ lp->tab = gbr_tab(bset, lp->row);
+ if (!lp->tab)
+ goto error;
+ lp->obj = NULL;
+ lp->neq = 0;
+
+ return lp;
+error:
+ delete_lp(lp);
+ return NULL;
+}
+
+static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
+{
+ lp->obj = row;
+}
+
+static int solve_lp(struct tab_lp *lp)
+{
+ enum isl_lp_result res;
+ unsigned flags = 0;
+
+ isl_int_set_si(lp->row->el[0], 0);
+ isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);
+ isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);
+ if (lp->neq)
+ flags = ISL_TAB_SAVE_DUAL;
+ res = isl_tab_min(lp->ctx, lp->tab, lp->row->el, lp->ctx->one,
+ &lp->opt, &lp->opt_denom, flags);
+ if (res != isl_lp_ok)
+ return -1;
+ return 0;
+}
+
+static void get_obj_val(struct tab_lp* lp, mpq_t *F)
+{
+ isl_int_neg(mpq_numref(*F), lp->opt);
+ isl_int_set(mpq_denref(*F), lp->opt_denom);
+}
+
+static void delete_lp(struct tab_lp *lp)
+{
+ if (!lp)
+ return;
+
+ isl_int_clear(lp->opt);
+ isl_int_clear(lp->opt_denom);
+ isl_vec_free(lp->row);
+ free(lp->stack);
+ isl_tab_free(lp->ctx, lp->tab);
+ isl_ctx_deref(lp->ctx);
+ free(lp);
+}
+
+static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
+{
+ lp->stack[lp->neq] = isl_tab_snap(lp->ctx, lp->tab);
+
+ isl_int_set_si(lp->row->el[0], 0);
+ isl_seq_cpy(lp->row->el + 1, row, lp->dim);
+ isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);
+
+ lp->tab = isl_tab_add_valid_eq(lp->ctx, lp->tab, lp->row->el);
+
+ return lp->neq++;
+}
+
+static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha)
+{
+ row += 2 * lp->n_ineq;
+ isl_int_neg(mpq_numref(*alpha), lp->tab->dual->el[1 + row]);
+ isl_int_set(mpq_denref(*alpha), lp->tab->dual->el[0]);
+}
+
+static void del_lp_row(struct tab_lp *lp)
+{
+ lp->neq--;
+ isl_tab_rollback(lp->ctx, lp->tab, lp->stack[lp->neq]);
+}
--- /dev/null
+#include <stdlib.h>
+#include "isl_basis_reduction.h"
+
+static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
+{
+ int i;
+
+ for (i = 0; i < n; ++i)
+ GBR_lp_get_alpha(lp, first + i, &alpha[i]);
+}
+
+/* This function implements the algorithm described in
+ * "An Implementation of the Generalized Basis Reduction Algorithm
+ * for Integer Programming" of Cook el al. to compute a reduced basis.
+ * We use \epsilon = 1/4.
+ *
+ * If options->gbr_only_first is set, the user is only interested
+ * in the first direction. In this case we stop the basis reduction when
+ * - the width in the first direction becomes smaller than 2
+ * or
+ * - we have moved forward all the way to the last direction
+ * and then back again all the way to the first direction.
+ */
+struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset)
+{
+ unsigned dim;
+ struct isl_mat *basis;
+ int unbounded;
+ int i;
+ GBR_LP *lp = NULL;
+ GBR_type F_old, alpha, F_new;
+ int row;
+ isl_int tmp;
+ struct isl_vec *b_tmp;
+ GBR_type *F = NULL;
+ GBR_type *alpha_buffer[2] = { NULL, NULL };
+ GBR_type *alpha_saved;
+ GBR_type F_saved;
+ int use_saved = 0;
+ isl_int mu[2];
+ GBR_type mu_F[2];
+ GBR_type two;
+
+ if (!bset)
+ return NULL;
+
+ dim = isl_basic_set_total_dim(bset);
+ basis = isl_mat_identity(bset->ctx, dim);
+ if (!basis)
+ return NULL;
+
+ if (dim == 1)
+ return basis;
+
+ isl_int_init(tmp);
+ isl_int_init(mu[0]);
+ isl_int_init(mu[1]);
+
+ GBR_init(alpha);
+ GBR_init(F_old);
+ GBR_init(F_new);
+ GBR_init(F_saved);
+ GBR_init(mu_F[0]);
+ GBR_init(mu_F[1]);
+ GBR_init(two);
+
+ b_tmp = isl_vec_alloc(bset->ctx, dim);
+ if (!b_tmp)
+ goto error;
+
+ F = isl_alloc_array(bset->ctx, GBR_type, dim);
+ alpha_buffer[0] = isl_alloc_array(bset->ctx, GBR_type, dim);
+ alpha_buffer[1] = isl_alloc_array(bset->ctx, GBR_type, dim);
+ alpha_saved = alpha_buffer[0];
+
+ if (!F || !alpha_buffer[0] || !alpha_buffer[1])
+ goto error;
+
+ for (i = 0; i < dim; ++i) {
+ GBR_init(F[i]);
+ GBR_init(alpha_buffer[0][i]);
+ GBR_init(alpha_buffer[1][i]);
+ }
+
+ GBR_set_ui(two, 2);
+
+ lp = GBR_lp_init(bset);
+ if (!lp)
+ goto error;
+
+ i = 0;
+
+ GBR_lp_set_obj(lp, basis->row[0], dim);
+ bset->ctx->stats->gbr_solved_lps++;
+ unbounded = GBR_lp_solve(lp);
+ isl_assert(bset->ctx, !unbounded, goto error);
+ GBR_lp_get_obj_val(lp, &F[0]);
+
+ do {
+ if (use_saved) {
+ row = GBR_lp_next_row(lp);
+ GBR_set(F_new, F_saved);
+ GBR_set(alpha, alpha_saved[i]);
+ } else {
+ row = GBR_lp_add_row(lp, basis->row[i], dim);
+ GBR_lp_set_obj(lp, basis->row[i+1], dim);
+ bset->ctx->stats->gbr_solved_lps++;
+ unbounded = GBR_lp_solve(lp);
+ isl_assert(bset->ctx, !unbounded, goto error);
+ GBR_lp_get_obj_val(lp, &F_new);
+
+ GBR_lp_get_alpha(lp, row, &alpha);
+
+ if (i > 0)
+ save_alpha(lp, row-i, i, alpha_saved);
+
+ GBR_lp_del_row(lp);
+ }
+ GBR_set(F[i+1], F_new);
+
+ GBR_floor(mu[0], alpha);
+ GBR_ceil(mu[1], alpha);
+
+ if (isl_int_eq(mu[0], mu[1]))
+ isl_int_set(tmp, mu[0]);
+ else {
+ int j;
+
+ for (j = 0; j <= 1; ++j) {
+ isl_int_set(tmp, mu[j]);
+ isl_seq_combine(b_tmp->el,
+ bset->ctx->one, basis->row[i+1],
+ tmp, basis->row[i], dim);
+ GBR_lp_set_obj(lp, b_tmp->el, dim);
+ bset->ctx->stats->gbr_solved_lps++;
+ unbounded = GBR_lp_solve(lp);
+ isl_assert(bset->ctx, !unbounded, goto error);
+ GBR_lp_get_obj_val(lp, &mu_F[j]);
+ if (i > 0)
+ save_alpha(lp, row-i, i, alpha_buffer[j]);
+ }
+
+ if (GBR_lt(mu_F[0], mu_F[1]))
+ j = 0;
+ else
+ j = 1;
+
+ isl_int_set(tmp, mu[j]);
+ GBR_set(F_new, mu_F[j]);
+ alpha_saved = alpha_buffer[j];
+ }
+ isl_seq_combine(basis->row[i+1],
+ bset->ctx->one, basis->row[i+1],
+ tmp, basis->row[i], dim);
+
+ GBR_set(F_old, F[i]);
+
+ use_saved = 0;
+ /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
+ GBR_set_ui(mu_F[0], 4);
+ GBR_mul(mu_F[0], mu_F[0], F_new);
+ GBR_set_ui(mu_F[1], 3);
+ GBR_mul(mu_F[1], mu_F[1], F_old);
+ if (GBR_lt(mu_F[0], mu_F[1])) {
+ basis = isl_mat_swap_rows(bset->ctx, basis, i, i + 1);
+ if (i > 0) {
+ use_saved = 1;
+ GBR_set(F_saved, F_new);
+ GBR_lp_del_row(lp);
+ --i;
+ } else {
+ GBR_set(F[0], F_new);
+ if (bset->ctx->gbr_only_first &&
+ GBR_lt(F[0], two))
+ break;
+ }
+ } else {
+ GBR_lp_add_row(lp, basis->row[i], dim);
+ ++i;
+ }
+ } while (i < dim-1);
+
+ if (0) {
+error:
+ isl_mat_free(bset->ctx, basis);
+ basis = NULL;
+ }
+
+ GBR_lp_delete(lp);
+
+ if (alpha_buffer[1])
+ for (i = 0; i < dim; ++i) {
+ GBR_clear(F[i]);
+ GBR_clear(alpha_buffer[0][i]);
+ GBR_clear(alpha_buffer[1][i]);
+ }
+ free(F);
+ free(alpha_buffer[0]);
+ free(alpha_buffer[1]);
+
+ isl_vec_free(b_tmp);
+
+ GBR_clear(alpha);
+ GBR_clear(F_old);
+ GBR_clear(F_new);
+ GBR_clear(F_saved);
+ GBR_clear(mu_F[0]);
+ GBR_clear(mu_F[1]);
+ GBR_clear(two);
+
+ isl_int_clear(tmp);
+ isl_int_clear(mu[0]);
+ isl_int_clear(mu[1]);
+
+ return basis;
+}
#include "isl_seq.h"
#include "isl_map_private.h"
#include "isl_equalities.h"
+#include "isl_tab.h"
+#include "isl_basis_reduction.h"
static struct isl_vec *empty_sample(struct isl_basic_set *bset)
{
return sample;
}
+/* Given a basic set "bset" and an affine function "f"/"denom",
+ * check if bset is bounded and non-empty and if so, return the minimal
+ * and maximal value attained by the affine function in "min" and "max".
+ * The minimal value is rounded up to the nearest integer, while the
+ * maximal value is rounded down.
+ * The return value indicates whether the set was empty or unbounded.
+ */
+static enum isl_lp_result basic_set_range(struct isl_basic_set *bset,
+ isl_int *f, isl_int denom, isl_int *min, isl_int *max)
+{
+ unsigned dim;
+ struct isl_tab *tab;
+ enum isl_lp_result res;
+
+ if (!bset)
+ return isl_lp_error;
+ if (isl_basic_set_fast_is_empty(bset))
+ return isl_lp_empty;
+
+ tab = isl_tab_from_basic_set(bset);
+ res = isl_tab_min(bset->ctx, tab, f, denom, min, NULL, 0);
+ if (res != isl_lp_ok) {
+ isl_tab_free(bset->ctx, tab);
+ return res;
+ }
+ dim = isl_basic_set_total_dim(bset);
+ isl_seq_neg(f, f, 1 + dim);
+ res = isl_tab_min(bset->ctx, tab, f, denom, max, NULL, 0);
+ isl_seq_neg(f, f, 1 + dim);
+ isl_int_neg(*max, *max);
+
+ isl_tab_free(bset->ctx, tab);
+ return res;
+}
+
+/* Perform a basis reduction on "bset" and return the inverse of
+ * the new basis, i.e., an affine mapping from the new coordinates to the old,
+ * in *T.
+ */
+static struct isl_basic_set *basic_set_reduced(struct isl_basic_set *bset,
+ struct isl_mat **T)
+{
+ struct isl_ctx *ctx;
+ unsigned gbr_only_first;
+
+ *T = NULL;
+ if (!bset)
+ return NULL;
+
+ ctx = bset->ctx;
+
+ gbr_only_first = ctx->gbr_only_first;
+ ctx->gbr_only_first = 1;
+ *T = isl_basic_set_reduced_basis(bset);
+ ctx->gbr_only_first = gbr_only_first;
+
+ *T = isl_mat_lin_to_aff(bset->ctx, *T);
+ *T = isl_mat_right_inverse(bset->ctx, *T);
+
+ bset = isl_basic_set_preimage(bset, isl_mat_copy(bset->ctx, *T));
+ if (!bset)
+ goto error;
+
+ return bset;
+error:
+ isl_mat_free(ctx, *T);
+ *T = NULL;
+ return NULL;
+}
+
+static struct isl_vec *sample_bounded(struct isl_basic_set *bset);
+
+/* Given a basic set "bset" whose first coordinate ranges between
+ * "min" and "max", step through all values from min to max, until
+ * the slice of bset with the first coordinate fixed to one of these
+ * values contains an integer point. If such a point is found, return it.
+ * If none of the slices contains any integer point, then bset itself
+ * doesn't contain any integer point and an empty sample is returned.
+ */
+static struct isl_vec *sample_scan(struct isl_basic_set *bset,
+ isl_int min, isl_int max)
+{
+ unsigned total;
+ struct isl_basic_set *slice = NULL;
+ struct isl_vec *sample = NULL;
+ isl_int tmp;
+
+ total = isl_basic_set_total_dim(bset);
+
+ isl_int_init(tmp);
+ for (isl_int_set(tmp, min); isl_int_le(tmp, max);
+ isl_int_add_ui(tmp, tmp, 1)) {
+ int k;
+
+ slice = isl_basic_set_copy(bset);
+ slice = isl_basic_set_cow(slice);
+ slice = isl_basic_set_extend_constraints(slice, 1, 0);
+ k = isl_basic_set_alloc_equality(slice);
+ if (k < 0)
+ goto error;
+ isl_int_set(slice->eq[k][0], tmp);
+ isl_int_set_si(slice->eq[k][1], -1);
+ isl_seq_clr(slice->eq[k] + 2, total - 1);
+ slice = isl_basic_set_simplify(slice);
+ sample = sample_bounded(slice);
+ slice = NULL;
+ if (!sample)
+ goto error;
+ if (sample->size > 0)
+ break;
+ isl_vec_free(sample);
+ sample = NULL;
+ }
+ if (!sample)
+ sample = empty_sample(bset);
+ else
+ isl_basic_set_free(bset);
+ isl_int_clear(tmp);
+ return sample;
+error:
+ isl_basic_set_free(bset);
+ isl_basic_set_free(slice);
+ isl_int_clear(tmp);
+ return NULL;
+}
+
+/* Given a basic set that is known to be bounded, find and return
+ * an integer point in the basic set, if there is any.
+ *
+ * After handling some trivial cases, we check the range of the
+ * first coordinate. If this coordinate can only attain one integer
+ * value, we are happy. Otherwise, we perform basis reduction and
+ * determine the new range.
+ *
+ * Then we step through all possible values in the range in sample_scan.
+ *
+ * If any basis reduction was performed, the sample value found, if any,
+ * is transformed back to the original space.
+ */
+static struct isl_vec *sample_bounded(struct isl_basic_set *bset)
+{
+ unsigned dim;
+ struct isl_ctx *ctx;
+ struct isl_vec *sample;
+ struct isl_vec *obj = NULL;
+ struct isl_mat *T = NULL;
+ isl_int min, max;
+ enum isl_lp_result res;
+
+ if (!bset)
+ return NULL;
+
+ if (isl_basic_set_fast_is_empty(bset))
+ return empty_sample(bset);
+
+ ctx = bset->ctx;
+ dim = isl_basic_set_total_dim(bset);
+ if (dim == 0)
+ return zero_sample(bset);
+ if (dim == 1)
+ return interval_sample(bset);
+ if (bset->n_eq > 0)
+ return sample_eq(bset, sample_bounded);
+
+ isl_int_init(min);
+ isl_int_init(max);
+ obj = isl_vec_alloc(bset->ctx, 1 + dim);
+ if (!obj)
+ goto error;
+ isl_seq_clr(obj->el, 1+ dim);
+ isl_int_set_si(obj->el[1], 1);
+
+ res = basic_set_range(bset, obj->el, bset->ctx->one, &min, &max);
+ if (res == isl_lp_error)
+ goto error;
+ isl_assert(bset->ctx, res != isl_lp_unbounded, goto error);
+ if (res == isl_lp_empty || isl_int_lt(max, min)) {
+ sample = empty_sample(bset);
+ goto out;
+ }
+
+ if (isl_int_ne(min, max)) {
+ bset = basic_set_reduced(bset, &T);
+ if (!bset)
+ goto error;
+
+ res = basic_set_range(bset, obj->el, bset->ctx->one, &min, &max);
+ if (res == isl_lp_error)
+ goto error;
+ isl_assert(bset->ctx, res != isl_lp_unbounded, goto error);
+ if (res == isl_lp_empty || isl_int_lt(max, min)) {
+ sample = empty_sample(bset);
+ goto out;
+ }
+ }
+
+ sample = sample_scan(bset, min, max);
+out:
+ if (T) {
+ if (!sample || sample->size == 0)
+ isl_mat_free(ctx, T);
+ else
+ sample = isl_mat_vec_product(ctx, T, sample);
+ }
+ isl_vec_free(obj);
+ isl_int_clear(min);
+ isl_int_clear(max);
+ return sample;
+error:
+ isl_mat_free(ctx, T);
+ isl_basic_set_free(bset);
+ isl_vec_free(obj);
+ isl_int_clear(min);
+ isl_int_clear(max);
+ return NULL;
+}
+
+/* Given a basic set "bset" and a value "sample" for the first coordinates
+ * of bset, plug in these values and drop the corresponding coordinates.
+ *
+ * We do this by computing the preimage of the transformation
+ *
+ * [ 1 0 ]
+ * x = [ s 0 ] x'
+ * [ 0 I ]
+ *
+ * where [1 s] is the sample value and I is the identity matrix of the
+ * appropriate dimension.
+ */
+static struct isl_basic_set *plug_in(struct isl_basic_set *bset,
+ struct isl_vec *sample)
+{
+ int i;
+ unsigned total;
+ struct isl_mat *T;
+
+ if (!bset || !sample)
+ goto error;
+
+ total = isl_basic_set_total_dim(bset);
+ T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
+ if (!T)
+ goto error;
+
+ for (i = 0; i < sample->size; ++i) {
+ isl_int_set(T->row[i][0], sample->el[i]);
+ isl_seq_clr(T->row[i] + 1, T->n_col - 1);
+ }
+ for (i = 0; i < T->n_col - 1; ++i) {
+ isl_seq_clr(T->row[sample->size + i], T->n_col);
+ isl_int_set_si(T->row[sample->size + i][1 + i], 1);
+ }
+ isl_vec_free(sample);
+
+ bset = isl_basic_set_preimage(bset, T);
+ return bset;
+error:
+ isl_basic_set_free(bset);
+ isl_vec_free(sample);
+ return NULL;
+}
+
+/* Given a basic set "bset", return any (possibly non-integer) point
+ * in the basic set.
+ */
+static struct isl_vec *rational_sample(struct isl_basic_set *bset)
+{
+ struct isl_tab *tab;
+ struct isl_vec *sample;
+
+ if (!bset)
+ return NULL;
+
+ tab = isl_tab_from_basic_set(bset);
+ sample = isl_tab_get_sample_value(bset->ctx, tab);
+ isl_tab_free(bset->ctx, tab);
+
+ isl_basic_set_free(bset);
+
+ return sample;
+}
+
+/* Given a rational vector, with the denominator in the first element
+ * of the vector, round up all coordinates.
+ */
+struct isl_vec *isl_vec_ceil(struct isl_vec *vec)
+{
+ int i;
+
+ vec = isl_vec_cow(vec);
+ if (!vec)
+ return NULL;
+
+ isl_seq_cdiv_q(vec->el + 1, vec->el + 1, vec->el[0], vec->size - 1);
+
+ isl_int_set_si(vec->el[0], 1);
+
+ return vec;
+}
+
+/* Given a linear cone "cone" and a rational point "vec",
+ * construct a polyhedron with shifted copies of the constraints in "cone",
+ * i.e., a polyhedron with "cone" as its recession cone, such that each
+ * point x in this polyhedron is such that the unit box positioned at x
+ * lies entirely inside the affine cone 'vec + cone'.
+ * Any rational point in this polyhedron may therefore be rounded up
+ * to yield an integer point that lies inside said affine cone.
+ *
+ * Denote the constraints of cone by "<a_i, x> >= 0" and the rational
+ * point "vec" by v/d.
+ * Let b_i = <a_i, v>. Then the affine cone 'vec + cone' is given
+ * by <a_i, x> - b/d >= 0.
+ * The polyhedron <a_i, x> - ceil{b/d} >= 0 is a subset of this affine cone.
+ * We prefer this polyhedron over the actual affine cone because it doesn't
+ * require a scaling of the constraints.
+ * If each of the vertices of the unit cube positioned at x lies inside
+ * this polyhedron, then the whole unit cube at x lies inside the affine cone.
+ * We therefore impose that x' = x + \sum e_i, for any selection of unit
+ * vectors lies inside the polyhedron, i.e.,
+ *
+ * <a_i, x'> - ceil{b/d} = <a_i, x> + sum a_i - ceil{b/d} >= 0
+ *
+ * The most stringent of these constraints is the one that selects
+ * all negative a_i, so the polyhedron we are looking for has constraints
+ *
+ * <a_i, x> + sum_{a_i < 0} a_i - ceil{b/d} >= 0
+ *
+ * Note that if cone were known to have only non-negative rays
+ * (which can be accomplished by a unimodular transformation),
+ * then we would only have to check the points x' = x + e_i
+ * and we only have to add the smallest negative a_i (if any)
+ * instead of the sum of all negative a_i.
+ */
+static struct isl_basic_set *shift_cone(struct isl_basic_set *cone,
+ struct isl_vec *vec)
+{
+ int i, j, k;
+ unsigned total;
+
+ struct isl_basic_set *shift = NULL;
+
+ if (!cone || !vec)
+ goto error;
+
+ isl_assert(cone->ctx, cone->n_eq == 0, goto error);
+
+ total = isl_basic_set_total_dim(cone);
+
+ shift = isl_basic_set_alloc_dim(isl_basic_set_get_dim(cone),
+ 0, 0, cone->n_ineq);
+
+ for (i = 0; i < cone->n_ineq; ++i) {
+ k = isl_basic_set_alloc_inequality(shift);
+ if (k < 0)
+ goto error;
+ isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
+ isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
+ &shift->ineq[k][0]);
+ isl_int_cdiv_q(shift->ineq[k][0],
+ shift->ineq[k][0], vec->el[0]);
+ isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
+ for (j = 0; j < total; ++j) {
+ if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
+ continue;
+ isl_int_add(shift->ineq[k][0],
+ shift->ineq[k][0], shift->ineq[k][1 + j]);
+ }
+ }
+
+ isl_basic_set_free(cone);
+ isl_vec_free(vec);
+
+ return isl_basic_set_finalize(shift);
+error:
+ isl_basic_set_free(shift);
+ isl_basic_set_free(cone);
+ isl_vec_free(vec);
+ return NULL;
+}
+
+/* Given a rational point vec in a (transformed) basic set,
+ * such that cone is the recession cone of the original basic set,
+ * "round up" the rational point to an integer point.
+ *
+ * We first check if the rational point just happens to be integer.
+ * If not, we transform the cone in the same way as the basic set,
+ * pick a point x in this cone shifted to the rational point such that
+ * the whole unit cube at x is also inside this affine cone.
+ * Then we simply round up the coordinates of x and return the
+ * resulting integer point.
+ */
+static struct isl_vec *round_up_in_cone(struct isl_vec *vec,
+ struct isl_basic_set *cone, struct isl_mat *U)
+{
+ unsigned total;
+
+ if (!vec || !cone || !U)
+ goto error;
+
+ isl_assert(vec->ctx, vec->size != 0, goto error);
+ if (isl_int_is_one(vec->el[0])) {
+ isl_mat_free(vec->ctx, U);
+ isl_basic_set_free(cone);
+ return vec;
+ }
+
+ total = isl_basic_set_total_dim(cone);
+ cone = isl_basic_set_preimage(cone, U);
+ cone = isl_basic_set_remove_dims(cone, 0, total - (vec->size - 1));
+
+ cone = shift_cone(cone, vec);
+
+ vec = rational_sample(cone);
+ vec = isl_vec_ceil(vec);
+ return vec;
+error:
+ isl_mat_free(vec ? vec->ctx : cone ? cone->ctx : NULL, U);
+ isl_vec_free(vec);
+ isl_basic_set_free(cone);
+ return NULL;
+}
+
+/* Concatenate two integer vectors, i.e., two vectors with denominator
+ * (stored in element 0) equal to 1.
+ */
+static struct isl_vec *vec_concat(struct isl_vec *vec1, struct isl_vec *vec2)
+{
+ struct isl_vec *vec;
+
+ if (!vec1 || !vec2)
+ goto error;
+ isl_assert(vec1->ctx, vec1->size > 0, goto error);
+ isl_assert(vec2->ctx, vec2->size > 0, goto error);
+ isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);
+ isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);
+
+ vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
+ if (!vec)
+ goto error;
+
+ isl_seq_cpy(vec->el, vec1->el, vec1->size);
+ isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
+
+ isl_vec_free(vec1);
+ isl_vec_free(vec2);
+
+ return vec;
+error:
+ isl_vec_free(vec1);
+ isl_vec_free(vec2);
+ return NULL;
+}
+
+/* Drop all constraints in bset that involve any of the dimensions
+ * first to first+n-1.
+ */
+static struct isl_basic_set *drop_constraints_involving
+ (struct isl_basic_set *bset, unsigned first, unsigned n)
+{
+ int i;
+
+ if (!bset)
+ return NULL;
+
+ bset = isl_basic_set_cow(bset);
+
+ for (i = bset->n_ineq - 1; i >= 0; --i) {
+ if (isl_seq_first_non_zero(bset->ineq[i] + 1 + first, n) == -1)
+ continue;
+ isl_basic_set_drop_inequality(bset, i);
+ }
+
+ return bset;
+}
+
+/* Give a basic set "bset" with recession cone "cone", compute and
+ * return an integer point in bset, if any.
+ *
+ * If the recession cone is full-dimensional, then we know that
+ * bset contains an infinite number of integer points and it is
+ * fairly easy to pick one of them.
+ * If the recession cone is not full-dimensional, then we first
+ * transform bset such that the bounded directions appear as
+ * the first dimensions of the transformed basic set.
+ * We do this by using a unimodular transformation that transforms
+ * the equalities in the recession cone to equalities on the first
+ * dimensions.
+ *
+ * The transformed set is then projected onto its bounded dimensions.
+ * Note that to compute this projection, we can simply drop all constraints
+ * involving any of the unbounded dimensions since these constraints
+ * cannot be combined to produce a constraint on the bounded dimensions.
+ * To see this, assume that there is such a combination of constraints
+ * that produces a constraint on the bounded dimensions. This means
+ * that some combination of the unbounded dimensions has both an upper
+ * bound and a lower bound in terms of the bounded dimensions, but then
+ * this combination would be a bounded direction too and would have been
+ * transformed into a bounded dimensions.
+ *
+ * We then compute a sample value in the bounded dimensions.
+ * If no such value can be found, then the original set did not contain
+ * any integer points and we are done.
+ * Otherwise, we plug in the value we found in the bounded dimensions,
+ * project out these bounded dimensions and end up with a set with
+ * a full-dimensional recession cone.
+ * A sample point in this set is computed by "rounding up" any
+ * rational point in the set.
+ *
+ * The sample points in the bounded and unbounded dimensions are
+ * then combined into a single sample point and transformed back
+ * to the original space.
+ */
+static struct isl_vec *sample_with_cone(struct isl_basic_set *bset,
+ struct isl_basic_set *cone)
+{
+ unsigned total;
+ unsigned cone_dim;
+ struct isl_mat *M, *U;
+ struct isl_vec *sample;
+ struct isl_vec *cone_sample;
+ struct isl_ctx *ctx;
+ struct isl_basic_set *bounded;
+
+ if (!bset || !cone)
+ goto error;
+
+ ctx = bset->ctx;
+ total = isl_basic_set_total_dim(cone);
+ cone_dim = total - cone->n_eq;
+
+ M = isl_mat_sub_alloc(bset->ctx, cone->eq, 0, cone->n_eq, 1, total);
+ M = isl_mat_left_hermite(bset->ctx, M, 0, &U, NULL);
+ if (!M)
+ goto error;
+ isl_mat_free(bset->ctx, M);
+
+ U = isl_mat_lin_to_aff(bset->ctx, U);
+ bset = isl_basic_set_preimage(bset, isl_mat_copy(bset->ctx, U));
+
+ bounded = isl_basic_set_copy(bset);
+ bounded = drop_constraints_involving(bounded, total - cone_dim, cone_dim);
+ bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
+ sample = sample_bounded(bounded);
+ if (!sample || sample->size == 0) {
+ isl_basic_set_free(bset);
+ isl_basic_set_free(cone);
+ isl_mat_free(ctx, U);
+ return sample;
+ }
+ bset = plug_in(bset, isl_vec_copy(sample));
+ cone_sample = rational_sample(bset);
+ cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(ctx, U));
+ sample = vec_concat(sample, cone_sample);
+ sample = isl_mat_vec_product(ctx, U, sample);
+ return sample;
+error:
+ isl_basic_set_free(cone);
+ isl_basic_set_free(bset);
+ return NULL;
+}
+
+/* Compute and return a sample point in bset using generalized basis
+ * reduction. We first check if the input set has a non-trivial
+ * recession cone. If so, we perform some extra preprocessing in
+ * sample_with_cone. Otherwise, we directly perform generalized basis
+ * reduction.
+ */
+static struct isl_vec *gbr_sample_no_lineality(struct isl_basic_set *bset)
+{
+ unsigned dim;
+ struct isl_basic_set *cone;
+
+ dim = isl_basic_set_total_dim(bset);
+
+ cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
+
+ if (cone->n_eq < dim)
+ return sample_with_cone(bset, cone);
+
+ isl_basic_set_free(cone);
+ return sample_bounded(bset);
+}
+
static struct isl_vec *sample_no_lineality(struct isl_basic_set *bset)
{
unsigned dim;
if (dim == 1)
return interval_sample(bset);
- return isl_pip_basic_set_sample(bset);
+ switch (bset->ctx->ilp_solver) {
+ case ISL_ILP_PIP:
+ return isl_pip_basic_set_sample(bset);
+ case ISL_ILP_GBR:
+ return gbr_sample_no_lineality(bset);
+ }
+ isl_assert(bset->ctx, 0, );
+ isl_basic_set_free(bset);
+ return NULL;
}
/* Compute an integer point in "bset" with a lineality space that