<commit_before>/*
 * Copyright 2008-2009 Katholieke Universiteit Leuven
 *
 * Use of this software is governed by the GNU LGPLv2.1 license
 *
 * Written by <NAME>, K.U.Leuven, Departement
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
 */

#include <isl_ctx_private.h>
#include <isl_map_private.h>
#include <isl/lp.h>
#include <isl/map.h>
#include <isl_mat_private.h>
#include <isl/set.h>
#include <isl/seq.h>
#include "isl_equalities.h"
#include "isl_tab.h"

static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set);

/* Return 1 if constraint c is redundant with respect to the constraints
 * in bmap.  If c is a lower [upper] bound in some variable and bmap
 * does not have a lower [upper] bound in that variable, then c cannot
 * be redundant and we do not need solve any lp.
 */
int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap,
	isl_int *c, isl_int *opt_n, isl_int *opt_d)
{
	enum isl_lp_result res;
	unsigned total;
	int i, j;

	if (!bmap)
		return -1;

	total = isl_basic_map_total_dim(*bmap);
	for (i = 0; i < total; ++i) {
		int sign;
		if (isl_int_is_zero(c[1+i]))
			continue;
		sign = isl_int_sgn(c[1+i]);
		for (j = 0; j < (*bmap)->n_ineq; ++j)
			if (sign == isl_int_sgn((*bmap)->ineq[j][1+i]))
				break;
		if (j == (*bmap)->n_ineq)
			break;
	}
	if (i < total)
		return 0;

	res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one,
					opt_n, opt_d, NULL);
	if (res == isl_lp_unbounded)
		return 0;
	if (res == isl_lp_error)
		return -1;
	if (res == isl_lp_empty) {
		*bmap = isl_basic_map_set_to_empty(*bmap);
		return 0;
	}
	return !isl_int_is_neg(*opt_n);
}

int isl_basic_set_constraint_is_redundant(struct isl_basic_set **bset,
	isl_int *c, isl_int *opt_n, isl_int *opt_d)
{
	return isl_basic_map_constraint_is_redundant(
			(struct isl_basic_map **)bset, c, opt_n, opt_d);
}

/* Remove redundant
 * constraints.  If the minimal value along the normal of a constraint
 * is the same if the constraint is removed, then the constraint is redundant.
 *
 * Alternatively, we could have intersected the basic map with the
 * corresponding equality and the checked if the dimension was that
 * of a facet.
 */
__isl_give isl_basic_map *isl_basic_map_remove_redundancies(
	__isl_take isl_basic_map *bmap)
{
	struct isl_tab *tab;

	if (!bmap)
		return NULL;

	bmap = isl_basic_map_gauss(bmap, NULL);
	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
		return bmap;
	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
		return bmap;
	if (bmap->n_ineq <= 1)
		return bmap;

	tab = isl_tab_from_basic_map(bmap);
	if (isl_tab_detect_implicit_equalities(tab) < 0)
		goto error;
	if (isl_tab_detect_redundant(tab) < 0)
		goto error;
	bmap = isl_basic_map_update_from_tab(bmap, tab);
	isl_tab_free(tab);
	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
	return bmap;
error:
	isl_tab_free(tab);
	isl_basic_map_free(bmap);
	return NULL;
}

__isl_give isl_basic_set *isl_basic_set_remove_redundancies(
	__isl_take isl_basic_set *bset)
{
	return (struct isl_basic_set *)
		isl_basic_map_remove_redundancies((struct isl_basic_map *)bset);
}

/* Check if the set set is bound in the direction of the affine
 * constraint c and if so, set the constant term such that the
 * resulting constraint is a bounding constraint for the set.
 */
static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
{
	int first;
	int j;
	isl_int opt;
	isl_int opt_denom;

	isl_int_init(opt);
	isl_int_init(opt_denom);
	first = 1;
	for (j = 0; j < set->n; ++j) {
		enum isl_lp_result res;

		if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
			continue;

		res = isl_basic_set_solve_lp(set->p[j],
				0, c, set->ctx->one, &opt, &opt_denom, NULL);
		if (res == isl_lp_unbounded)
			break;
		if (res == isl_lp_error)
			goto error;
		if (res == isl_lp_empty) {
			set->p[j] = isl_basic_set_set_to_empty(set->p[j]);
			if (!set->p[j])
				goto error;
			continue;
		}
		if (first || isl_int_is_neg(opt)) {
			if (!isl_int_is_one(opt_denom))
				isl_seq_scale(c, c, opt_denom, len);
			isl_int_sub(c[0], c[0], opt);
		}
		first = 0;
	}
	isl_int_clear(opt);
	isl_int_clear(opt_denom);
	return j >= set->n;
error:
	isl_int_clear(opt);
	isl_int_clear(opt_denom);
	return -1;
}

__isl_give isl_basic_map *isl_basic_map_set_rational(
	__isl_take isl_basic_set *bmap)
{
	if (!bmap)
		return NULL;

	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
		return bmap;

	bmap = isl_basic_map_cow(bmap);
	if (!bmap)
		return NULL;

	ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);

	return isl_basic_map_finalize(bmap);
}

__isl_give isl_basic_set *isl_basic_set_set_rational(
	__isl_take isl_basic_set *bset)
{
	return isl_basic_map_set_rational(bset);
}

static struct isl_set *isl_set_set_rational(struct isl_set *set)
{
	int i;

	set = isl_set_cow(set);
	if (!set)
		return NULL;
	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_set_rational(set->p[i]);
		if (!set->p[i])
			goto error;
	}
	return set;
error:
	isl_set_free(set);
	return NULL;
}

static struct isl_basic_set *isl_basic_set_add_equality(
	struct isl_basic_set *bset, isl_int *c)
{
	int i;
	unsigned dim;

	if (!bset)
		return NULL;

	if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
		return bset;

	isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
	isl_assert(bset->ctx, bset->n_div == 0, goto error);
	dim = isl_basic_set_n_dim(bset);
	bset = isl_basic_set_cow(bset);
	bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
	i = isl_basic_set_alloc_equality(bset);
	if (i < 0)
		goto error;
	isl_seq_cpy(bset->eq[i], c, 1 + dim);
	return bset;
error:
	isl_basic_set_free(bset);
	return NULL;
}

static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c)
{
	int i;

	set = isl_set_cow(set);
	if (!set)
		return NULL;
	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_add_equality(set->p[i], c);
		if (!set->p[i])
			goto error;
	}
	return set;
error:
	isl_set_free(set);
	return NULL;
}

/* Given a union of basic sets, construct the constraints for wrapping
 * a facet around one of its ridges.
 * In particular, if each of n the d-dimensional basic sets i in "set"
 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
 * and is defined by the constraints
 *				    [ 1 ]
 *				A_i [ x ]  >= 0
 *
 * then the resulting set is of dimension n*(1+d) and has as constraints
 *
 *				    [ a_i ]
 *				A_i [ x_i ] >= 0
 *
 *				      a_i   >= 0
 *
 *			\sum_i x_{i,1} = 1
 */
static struct isl_basic_set *wrap_constraints(struct isl_set *set)
{
	struct isl_basic_set *lp;
	unsigned n_eq;
	unsigned n_ineq;
	int i, j, k;
	unsigned dim, lp_dim;

	if (!set)
		return NULL;

	dim = 1 + isl_set_n_dim(set);
	n_eq = 1;
	n_ineq = set->n;
	for (i = 0; i < set->n; ++i) {
		n_eq += set->p[i]->n_eq;
		n_ineq += set->p[i]->n_ineq;
	}
	lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
	lp = isl_basic_set_set_rational(lp);
	if (!lp)
		return NULL;
	lp_dim = isl_basic_set_n_dim(lp);
	k = isl_basic_set_alloc_equality(lp);
	isl_int_set_si(lp->eq[k][0], -1);
	for (i = 0; i < set->n; ++i) {
		isl_int_set_si(lp->eq[k][1+dim*i], 0);
		isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
		isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
	}
	for (i = 0; i < set->n; ++i) {
		k = isl_basic_set_alloc_inequality(lp);
		isl_seq_clr(lp->ineq[k], 1+lp_dim);
		isl_int_set_si(lp->ineq[k][1+dim*i], 1);

		for (j = 0; j < set->p[i]->n_eq; ++j) {
			k = isl_basic_set_alloc_equality(lp);
			isl_seq_clr(lp->eq[k], 1+dim*i);
			isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
			isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
		}

		for (j = 0; j < set->p[i]->n_ineq; ++j) {
			k = isl_basic_set_alloc_inequality(lp);
			isl_seq_clr(lp->ineq[k], 1+dim*i);
			isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
			isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
		}
	}
	return lp;
}

/* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
 * of that facet, compute the other facet of the convex hull that contains
 * the ridge.
 *
 * We first transform the set such that the facet constraint becomes
 *
 *			x_1 >= 0
 *
 * I.e., the facet lies in
 *
 *			x_1 = 0
 *
 * and on that facet, the constraint that defines the ridge is
 *
 *			x_2 >= 0
 *
 * (This transformation is not strictly needed, all that is needed is
 * that the ridge contains the origin.)
 *
 * Since the ridge contains the origin, the cone of the convex hull
 * will be of the form
 *
 *			x_1 >= 0
 *			x_2 >= a x_1
 *
 * with this second constraint defining the new facet.
 * The constant a is obtained by settting x_1 in the cone of the
 * convex hull to 1 and minimizing x_2.
 * Now, each element in the cone of the convex hull is the sum
 * of elements in the cones of the basic sets.
 * If a_i is the dilation factor of basic set i, then the problem
 * we need to solve is
 *
 *			min \sum_i x_{i,2}
 *			st
 *				\sum_i x_{i,1} = 1
 *				    a_i   >= 0
 *				  [ a_i ]
 *				A [ x_i ] >= 0
 *
 * with
 *				    [  1  ]
 *				A_i [ x_i ] >= 0
 *
 * the constraints of each (transformed) basic set.
 * If a = n/d, then the constraint defining the new facet (in the transformed
 * space) is
 *
 *			-n x_1 + d x_2 >= 0
 *
 * In the original space, we need to take the same combination of the
 * corresponding constraints "facet" and "ridge".
 *
 * If a = -infty = "-1/0", then we just return the original facet constraint.
 * This means that the facet is unbounded, but has a bounded intersection
 * with the union of sets.
 */
isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
	isl_int *facet, isl_int *ridge)
{
	int i;
	isl_ctx *ctx;
	struct isl_mat *T = NULL;
	struct isl_basic_set *lp = NULL;
	struct isl_vec *obj;
	enum isl_lp_result res;
	isl_int num, den;
	unsigned dim;

	if (!set)
		return NULL;
	ctx = set->ctx;
	set = isl_set_copy(set);
	set = isl_set_set_rational(set);

	dim = 1 + isl_set_n_dim(set);
	T = isl_mat_alloc(ctx, 3, dim);
	if (!T)
		goto error;
	isl_int_set_si(T->row[0][0], 1);
	isl_seq_clr(T->row[0]+1, dim - 1);
	isl_seq_cpy(T->row[1], facet, dim);
	isl_seq_cpy(T->row[2], ridge, dim);
	T = isl_mat_right_inverse(T);
	set = isl_set_preimage(set, T);
	T = NULL;
	if (!set)
		goto error;
	lp = wrap_constraints(set);
	obj = isl_vec_alloc(ctx, 1 + dim*set->n);
	if (!obj)
		goto error;
	isl_int_set_si(obj->block.data[0], 0);
	for (i = 0; i < set->n; ++i) {
		isl_seq_clr(obj->block.data + 1 + dim*i, 2);
		isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
		isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
	}
	isl_int_init(num);
	isl_int_init(den);
	res = isl_basic_set_solve_lp(lp, 0,
			    obj->block.data, ctx->one, &num, &den, NULL);
	if (res == isl_lp_ok) {
		isl_int_neg(num, num);
		isl_seq_combine(facet, num, facet, den, ridge, dim);
		isl_seq_normalize(ctx, facet, dim);
	}
	isl_int_clear(num);
	isl_int_clear(den);
	isl_vec_free(obj);
	isl_basic_set_free(lp);
	isl_set_free(set);
	if (res == isl_lp_error)
		return NULL;
	isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded, 
		   return NULL);
	return facet;
error:
	isl_basic_set_free(lp);
	isl_mat_free(T);
	isl_set_free(set);
	return NULL;
}

/* Compute the constraint of a facet of "set".
 *
 * We first compute the intersection with a bounding constraint
 * that is orthogonal to one of the coordinate axes.
 * If the affine hull of this intersection has only one equality,
 * we have found a facet.
 * Otherwise, we wrap the current bounding constraint around
 * one of the equalities of the face (one that is not equal to
 * the current bounding constraint).
 * This process continues until we have found a facet.
 * The dimension of the intersection increases by at least
 * one on each iteration, so termination is guaranteed.
 */
static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
{
	struct isl_set *slice = NULL;
	struct isl_basic_set *face = NULL;
	int i;
	unsigned dim = isl_set_n_dim(set);
	int is_bound;
	isl_mat *bounds;

	isl_assert(set->ctx, set->n > 0, goto error);
	bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
	if (!bounds)
		return NULL;

	isl_seq_clr(bounds->row[0], dim);
	isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
	is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
	if (is_bound < 0)
		goto error;
	isl_assert(set->ctx, is_bound, goto error);
	isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
	bounds->n_row = 1;

	for (;;) {
		slice = isl_set_copy(set);
		slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
		face = isl_set_affine_hull(slice);
		if (!face)
			goto error;
		if (face->n_eq == 1) {
			isl_basic_set_free(face);
			break;
		}
		for (i = 0; i < face->n_eq; ++i)
			if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
			    !isl_seq_is_neg(bounds->row[0],
						face->eq[i], 1 + dim))
				break;
		isl_assert(set->ctx, i < face->n_eq, goto error);
		if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
			goto error;
		isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
		isl_basic_set_free(face);
	}

	return bounds;
error:
	isl_basic_set_free(face);
	isl_mat_free(bounds);
	return NULL;
}

/* Given the bounding constraint "c" of a facet of the convex hull of "set",
 * compute a hyperplane description of the facet, i.e., compute the facets
 * of the facet.
 *
 * We compute an affine transformation that transforms the constraint
 *
 *			  [ 1 ]
 *			c [ x ] = 0
 *
 * to the constraint
 *
 *			   z_1  = 0
 *
 * by computing the right inverse U of a matrix that starts with the rows
 *
 *			[ 1 0 ]
 *			[  c  ]
 *
 * Then
 *			[ 1 ]     [ 1 ]
 *			[ x ] = U [ z ]
 * and
 *			[ 1 ]     [ 1 ]
 *			[ z ] = Q [ x ]
 *
 * with Q = U^{-1}
 * Since z_1 is zero, we can drop this variable as well as the corresponding
 * column of U to obtain
 *
 *			[ 1 ]      [ 1  ]
 *			[ x ] = U' [ z' ]
 * and
 *			[ 1  ]      [ 1 ]
 *			[ z' ] = Q' [ x ]
 *
 * with Q' equal to Q, but without the corresponding row.
 * After computing the facets of the facet in the z' space,
 * we convert them back to the x space through Q.
 */
static struct isl_basic_set *compute_facet(struct isl_set *set, isl_int *c)
{
	struct isl_mat *m, *U, *Q;
	struct isl_basic_set *facet = NULL;
	struct isl_ctx *ctx;
	unsigned dim;

	ctx = set->ctx;
	set = isl_set_copy(set);
	dim = isl_set_n_dim(set);
	m = isl_mat_alloc(set->ctx, 2, 1 + dim);
	if (!m)
		goto error;
	isl_int_set_si(m->row[0][0], 1);
	isl_seq_clr(m->row[0]+1, dim);
	isl_seq_cpy(m->row[1], c, 1+dim);
	U = isl_mat_right_inverse(m);
	Q = isl_mat_right_inverse(isl_mat_copy(U));
	U = isl_mat_drop_cols(U, 1, 1);
	Q = isl_mat_drop_rows(Q, 1, 1);
	set = isl_set_preimage(set, U);
	facet = uset_convex_hull_wrap_bounded(set);
	facet = isl_basic_set_preimage(facet, Q);
	if (facet)
		isl_assert(ctx, facet->n_eq == 0, goto error);
	return facet;
error:
	isl_basic_set_free(facet);
	isl_set_free(set);
	return NULL;
}

/* Given an initial facet constraint, compute the remaining facets.
 * We do this by running through all facets found so far and computing
 * the adjacent facets through wrapping, adding those facets that we
 * hadn't already found before.
 *
 * For each facet we have found so far, we first compute its facets
 * in the resulting convex hull.  That is, we compute the ridges
 * of the resulting convex hull contained in the facet.
 * We also compute the corresponding facet in the current approximation
 * of the convex hull.  There is no need to wrap around the ridges
 * in this facet since that would result in a facet that is already
 * present in the current approximation.
 *
 * This function can still be significantly optimized by checking which of
 * the facets of the basic sets are also facets of the convex hull and
 * using all the facets so far to help in constructing the facets of the
 * facets
 * and/or
 * using the technique in section "3.1 Ridge Generation" of
 * "Extended Convex Hull" by Fukuda et al.
 */
static struct isl_basic_set *extend(struct isl_basic_set *hull,
	struct isl_set *set)
{
	int i, j, f;
	int k;
	struct isl_basic_set *facet = NULL;
	struct isl_basic_set *hull_facet = NULL;
	unsigned dim;

	if (!hull)
		return NULL;

	isl_assert(set->ctx, set->n > 0, goto error);

	dim = isl_set_n_dim(set);

	for (i = 0; i < hull->n_ineq; ++i) {
		facet = compute_facet(set, hull->ineq[i]);
		facet = isl_basic_set_add_equality(facet, hull->ineq[i]);
		facet = isl_basic_set_gauss(facet, NULL);
		facet = isl_basic_set_normalize_constraints(facet);
		hull_facet = isl_basic_set_copy(hull);
		hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]);
		hull_facet = isl_basic_set_gauss(hull_facet, NULL);
		hull_facet = isl_basic_set_normalize_constraints(hull_facet);
		if (!facet || !hull_facet)
			goto error;
		hull = isl_basic_set_cow(hull);
		hull = isl_basic_set_extend_dim(hull,
			isl_dim_copy(hull->dim), 0, 0, facet->n_ineq);
		if (!hull)
			goto error;
		for (j = 0; j < facet->n_ineq; ++j) {
			for (f = 0; f < hull_facet->n_ineq; ++f)
				if (isl_seq_eq(facet->ineq[j],
						hull_facet->ineq[f], 1 + dim))
					break;
			if (f < hull_facet->n_ineq)
				continue;
			k = isl_basic_set_alloc_inequality(hull);
			if (k < 0)
				goto error;
			isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
			if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
				goto error;
		}
		isl_basic_set_free(hull_facet);
		isl_basic_set_free(facet);
	}
	hull = isl_basic_set_simplify(hull);
	hull = isl_basic_set_finalize(hull);
	return hull;
error:
	isl_basic_set_free(hull_facet);
	isl_basic_set_free(facet);
	isl_basic_set_free(hull);
	return NULL;
}

/* Special case for computing the convex hull of a one dimensional set.
 * We simply collect the lower and upper bounds of each basic set
 * and the biggest of those.
 */
static struct isl_basic_set *convex_hull_1d(struct isl_set *set)
{
	struct isl_mat *c = NULL;
	isl_int *lower = NULL;
	isl_int *upper = NULL;
	int i, j, k;
	isl_int a, b;
	struct isl_basic_set *hull;

	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_simplify(set->p[i]);
		if (!set->p[i])
			goto error;
	}
	set = isl_set_remove_empty_parts(set);
	if (!set)
		goto error;
	isl_assert(set->ctx, set->n > 0, goto error);
	c = isl_mat_alloc(set->ctx, 2, 2);
	if (!c)
		goto error;

	if (set->p[0]->n_eq > 0) {
		isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
		lower = c->row[0];
		upper = c->row[1];
		if (isl_int_is_pos(set->p[0]->eq[0][1])) {
			isl_seq_cpy(lower, set->p[0]->eq[0], 2);
			isl_seq_neg(upper, set->p[0]->eq[0], 2);
		} else {
			isl_seq_neg(lower, set->p[0]->eq[0], 2);
			isl_seq_cpy(upper, set->p[0]->eq[0], 2);
		}
	} else {
		for (j = 0; j < set->p[0]->n_ineq; ++j) {
			if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
				lower = c->row[0];
				isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
			} else {
				upper = c->row[1];
				isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
			}
		}
	}

	isl_int_init(a);
	isl_int_init(b);
	for (i = 0; i < set->n; ++i) {
		struct isl_basic_set *bset = set->p[i];
		int has_lower = 0;
		int has_upper = 0;

		for (j = 0; j < bset->n_eq; ++j) {
			has_lower = 1;
			has_upper = 1;
			if (lower) {
				isl_int_mul(a, lower[0], bset->eq[j][1]);
				isl_int_mul(b, lower[1], bset->eq[j][0]);
				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
					isl_seq_cpy(lower, bset->eq[j], 2);
				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
					isl_seq_neg(lower, bset->eq[j], 2);
			}
			if (upper) {
				isl_int_mul(a, upper[0], bset->eq[j][1]);
				isl_int_mul(b, upper[1], bset->eq[j][0]);
				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
					isl_seq_neg(upper, bset->eq[j], 2);
				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
					isl_seq_cpy(upper, bset->eq[j], 2);
			}
		}
		for (j = 0; j < bset->n_ineq; ++j) {
			if (isl_int_is_pos(bset->ineq[j][1]))
				has_lower = 1;
			if (isl_int_is_neg(bset->ineq[j][1]))
				has_upper = 1;
			if (lower && isl_int_is_pos(bset->ineq[j][1])) {
				isl_int_mul(a, lower[0], bset->ineq[j][1]);
				isl_int_mul(b, lower[1], bset->ineq[j][0]);
				if (isl_int_lt(a, b))
					isl_seq_cpy(lower, bset->ineq[j], 2);
			}
			if (upper && isl_int_is_neg(bset->ineq[j][1])) {
				isl_int_mul(a, upper[0], bset->ineq[j][1]);
				isl_int_mul(b, upper[1], bset->ineq[j][0]);
				if (isl_int_gt(a, b))
					isl_seq_cpy(upper, bset->ineq[j], 2);
			}
		}
		if (!has_lower)
			lower = NULL;
		if (!has_upper)
			upper = NULL;
	}
	isl_int_clear(a);
	isl_int_clear(b);

	hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
	hull = isl_basic_set_set_rational(hull);
	if (!hull)
		goto error;
	if (lower) {
		k = isl_basic_set_alloc_inequality(hull);
		isl_seq_cpy(hull->ineq[k], lower, 2);
	}
	if (upper) {
		k = isl_basic_set_alloc_inequality(hull);
		isl_seq_cpy(hull->ineq[k], upper, 2);
	}
	hull = isl_basic_set_finalize(hull);
	isl_set_free(set);
	isl_mat_free(c);
	return hull;
error:
	isl_set_free(set);
	isl_mat_free(c);
	return NULL;
}

/* Project out final n dimensions using Fourier-Motzkin */
static struct isl_set *set_project_out(struct isl_ctx *ctx,
	struct isl_set *set, unsigned n)
{
	return isl_set_remove_dims(set, isl_dim_set, isl_set_n_dim(set) - n, n);
}

static struct isl_basic_set *convex_hull_0d(struct isl_set *set)
{
	struct isl_basic_set *convex_hull;

	if (!set)
		return NULL;

	if (isl_set_is_empty(set))
		convex_hull = isl_basic_set_empty(isl_dim_copy(set->dim));
	else
		convex_hull = isl_basic_set_universe(isl_dim_copy(set->dim));
	isl_set_free(set);
	return convex_hull;
}

/* Compute the convex hull of a pair of basic sets without any parameters or
 * integer divisions using Fourier-Motzkin elimination.
 * The convex hull is the set of all points that can be written as
 * the sum of points from both basic sets (in homogeneous coordinates).
 * We set up the constraints in a space with dimensions for each of
 * the three sets and then project out the dimensions corresponding
 * to the two original basic sets, retaining only those corresponding
 * to the convex hull.
 */
static struct isl_basic_set *convex_hull_pair_elim(struct isl_basic_set *bset1,
	struct isl_basic_set *bset2)
{
	int i, j, k;<commit_msg>isl_convex_hull.c: remove unused set_project_out<commit_after>/*
 * Copyright 2008-2009 Katholieke Universiteit Leuven
 *
 * Use of this software is governed by the GNU LGPLv2.1 license
 *
 * Written by <NAME>, K.U.Leuven, Departement
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
 */

#include <isl_ctx_private.h>
#include <isl_map_private.h>
#include <isl/lp.h>
#include <isl/map.h>
#include <isl_mat_private.h>
#include <isl/set.h>
#include <isl/seq.h>
#include "isl_equalities.h"
#include "isl_tab.h"

static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set);

/* Return 1 if constraint c is redundant with respect to the constraints
 * in bmap.  If c is a lower [upper] bound in some variable and bmap
 * does not have a lower [upper] bound in that variable, then c cannot
 * be redundant and we do not need solve any lp.
 */
int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap,
	isl_int *c, isl_int *opt_n, isl_int *opt_d)
{
	enum isl_lp_result res;
	unsigned total;
	int i, j;

	if (!bmap)
		return -1;

	total = isl_basic_map_total_dim(*bmap);
	for (i = 0; i < total; ++i) {
		int sign;
		if (isl_int_is_zero(c[1+i]))
			continue;
		sign = isl_int_sgn(c[1+i]);
		for (j = 0; j < (*bmap)->n_ineq; ++j)
			if (sign == isl_int_sgn((*bmap)->ineq[j][1+i]))
				break;
		if (j == (*bmap)->n_ineq)
			break;
	}
	if (i < total)
		return 0;

	res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one,
					opt_n, opt_d, NULL);
	if (res == isl_lp_unbounded)
		return 0;
	if (res == isl_lp_error)
		return -1;
	if (res == isl_lp_empty) {
		*bmap = isl_basic_map_set_to_empty(*bmap);
		return 0;
	}
	return !isl_int_is_neg(*opt_n);
}

int isl_basic_set_constraint_is_redundant(struct isl_basic_set **bset,
	isl_int *c, isl_int *opt_n, isl_int *opt_d)
{
	return isl_basic_map_constraint_is_redundant(
			(struct isl_basic_map **)bset, c, opt_n, opt_d);
}

/* Remove redundant
 * constraints.  If the minimal value along the normal of a constraint
 * is the same if the constraint is removed, then the constraint is redundant.
 *
 * Alternatively, we could have intersected the basic map with the
 * corresponding equality and the checked if the dimension was that
 * of a facet.
 */
__isl_give isl_basic_map *isl_basic_map_remove_redundancies(
	__isl_take isl_basic_map *bmap)
{
	struct isl_tab *tab;

	if (!bmap)
		return NULL;

	bmap = isl_basic_map_gauss(bmap, NULL);
	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
		return bmap;
	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
		return bmap;
	if (bmap->n_ineq <= 1)
		return bmap;

	tab = isl_tab_from_basic_map(bmap);
	if (isl_tab_detect_implicit_equalities(tab) < 0)
		goto error;
	if (isl_tab_detect_redundant(tab) < 0)
		goto error;
	bmap = isl_basic_map_update_from_tab(bmap, tab);
	isl_tab_free(tab);
	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
	return bmap;
error:
	isl_tab_free(tab);
	isl_basic_map_free(bmap);
	return NULL;
}

__isl_give isl_basic_set *isl_basic_set_remove_redundancies(
	__isl_take isl_basic_set *bset)
{
	return (struct isl_basic_set *)
		isl_basic_map_remove_redundancies((struct isl_basic_map *)bset);
}

/* Check if the set set is bound in the direction of the affine
 * constraint c and if so, set the constant term such that the
 * resulting constraint is a bounding constraint for the set.
 */
static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
{
	int first;
	int j;
	isl_int opt;
	isl_int opt_denom;

	isl_int_init(opt);
	isl_int_init(opt_denom);
	first = 1;
	for (j = 0; j < set->n; ++j) {
		enum isl_lp_result res;

		if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
			continue;

		res = isl_basic_set_solve_lp(set->p[j],
				0, c, set->ctx->one, &opt, &opt_denom, NULL);
		if (res == isl_lp_unbounded)
			break;
		if (res == isl_lp_error)
			goto error;
		if (res == isl_lp_empty) {
			set->p[j] = isl_basic_set_set_to_empty(set->p[j]);
			if (!set->p[j])
				goto error;
			continue;
		}
		if (first || isl_int_is_neg(opt)) {
			if (!isl_int_is_one(opt_denom))
				isl_seq_scale(c, c, opt_denom, len);
			isl_int_sub(c[0], c[0], opt);
		}
		first = 0;
	}
	isl_int_clear(opt);
	isl_int_clear(opt_denom);
	return j >= set->n;
error:
	isl_int_clear(opt);
	isl_int_clear(opt_denom);
	return -1;
}

__isl_give isl_basic_map *isl_basic_map_set_rational(
	__isl_take isl_basic_set *bmap)
{
	if (!bmap)
		return NULL;

	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
		return bmap;

	bmap = isl_basic_map_cow(bmap);
	if (!bmap)
		return NULL;

	ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL);

	return isl_basic_map_finalize(bmap);
}

__isl_give isl_basic_set *isl_basic_set_set_rational(
	__isl_take isl_basic_set *bset)
{
	return isl_basic_map_set_rational(bset);
}

static struct isl_set *isl_set_set_rational(struct isl_set *set)
{
	int i;

	set = isl_set_cow(set);
	if (!set)
		return NULL;
	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_set_rational(set->p[i]);
		if (!set->p[i])
			goto error;
	}
	return set;
error:
	isl_set_free(set);
	return NULL;
}

static struct isl_basic_set *isl_basic_set_add_equality(
	struct isl_basic_set *bset, isl_int *c)
{
	int i;
	unsigned dim;

	if (!bset)
		return NULL;

	if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
		return bset;

	isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
	isl_assert(bset->ctx, bset->n_div == 0, goto error);
	dim = isl_basic_set_n_dim(bset);
	bset = isl_basic_set_cow(bset);
	bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
	i = isl_basic_set_alloc_equality(bset);
	if (i < 0)
		goto error;
	isl_seq_cpy(bset->eq[i], c, 1 + dim);
	return bset;
error:
	isl_basic_set_free(bset);
	return NULL;
}

static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c)
{
	int i;

	set = isl_set_cow(set);
	if (!set)
		return NULL;
	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_add_equality(set->p[i], c);
		if (!set->p[i])
			goto error;
	}
	return set;
error:
	isl_set_free(set);
	return NULL;
}

/* Given a union of basic sets, construct the constraints for wrapping
 * a facet around one of its ridges.
 * In particular, if each of n the d-dimensional basic sets i in "set"
 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
 * and is defined by the constraints
 *				    [ 1 ]
 *				A_i [ x ]  >= 0
 *
 * then the resulting set is of dimension n*(1+d) and has as constraints
 *
 *				    [ a_i ]
 *				A_i [ x_i ] >= 0
 *
 *				      a_i   >= 0
 *
 *			\sum_i x_{i,1} = 1
 */
static struct isl_basic_set *wrap_constraints(struct isl_set *set)
{
	struct isl_basic_set *lp;
	unsigned n_eq;
	unsigned n_ineq;
	int i, j, k;
	unsigned dim, lp_dim;

	if (!set)
		return NULL;

	dim = 1 + isl_set_n_dim(set);
	n_eq = 1;
	n_ineq = set->n;
	for (i = 0; i < set->n; ++i) {
		n_eq += set->p[i]->n_eq;
		n_ineq += set->p[i]->n_ineq;
	}
	lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
	lp = isl_basic_set_set_rational(lp);
	if (!lp)
		return NULL;
	lp_dim = isl_basic_set_n_dim(lp);
	k = isl_basic_set_alloc_equality(lp);
	isl_int_set_si(lp->eq[k][0], -1);
	for (i = 0; i < set->n; ++i) {
		isl_int_set_si(lp->eq[k][1+dim*i], 0);
		isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
		isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
	}
	for (i = 0; i < set->n; ++i) {
		k = isl_basic_set_alloc_inequality(lp);
		isl_seq_clr(lp->ineq[k], 1+lp_dim);
		isl_int_set_si(lp->ineq[k][1+dim*i], 1);

		for (j = 0; j < set->p[i]->n_eq; ++j) {
			k = isl_basic_set_alloc_equality(lp);
			isl_seq_clr(lp->eq[k], 1+dim*i);
			isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
			isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
		}

		for (j = 0; j < set->p[i]->n_ineq; ++j) {
			k = isl_basic_set_alloc_inequality(lp);
			isl_seq_clr(lp->ineq[k], 1+dim*i);
			isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
			isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
		}
	}
	return lp;
}

/* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
 * of that facet, compute the other facet of the convex hull that contains
 * the ridge.
 *
 * We first transform the set such that the facet constraint becomes
 *
 *			x_1 >= 0
 *
 * I.e., the facet lies in
 *
 *			x_1 = 0
 *
 * and on that facet, the constraint that defines the ridge is
 *
 *			x_2 >= 0
 *
 * (This transformation is not strictly needed, all that is needed is
 * that the ridge contains the origin.)
 *
 * Since the ridge contains the origin, the cone of the convex hull
 * will be of the form
 *
 *			x_1 >= 0
 *			x_2 >= a x_1
 *
 * with this second constraint defining the new facet.
 * The constant a is obtained by settting x_1 in the cone of the
 * convex hull to 1 and minimizing x_2.
 * Now, each element in the cone of the convex hull is the sum
 * of elements in the cones of the basic sets.
 * If a_i is the dilation factor of basic set i, then the problem
 * we need to solve is
 *
 *			min \sum_i x_{i,2}
 *			st
 *				\sum_i x_{i,1} = 1
 *				    a_i   >= 0
 *				  [ a_i ]
 *				A [ x_i ] >= 0
 *
 * with
 *				    [  1  ]
 *				A_i [ x_i ] >= 0
 *
 * the constraints of each (transformed) basic set.
 * If a = n/d, then the constraint defining the new facet (in the transformed
 * space) is
 *
 *			-n x_1 + d x_2 >= 0
 *
 * In the original space, we need to take the same combination of the
 * corresponding constraints "facet" and "ridge".
 *
 * If a = -infty = "-1/0", then we just return the original facet constraint.
 * This means that the facet is unbounded, but has a bounded intersection
 * with the union of sets.
 */
isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
	isl_int *facet, isl_int *ridge)
{
	int i;
	isl_ctx *ctx;
	struct isl_mat *T = NULL;
	struct isl_basic_set *lp = NULL;
	struct isl_vec *obj;
	enum isl_lp_result res;
	isl_int num, den;
	unsigned dim;

	if (!set)
		return NULL;
	ctx = set->ctx;
	set = isl_set_copy(set);
	set = isl_set_set_rational(set);

	dim = 1 + isl_set_n_dim(set);
	T = isl_mat_alloc(ctx, 3, dim);
	if (!T)
		goto error;
	isl_int_set_si(T->row[0][0], 1);
	isl_seq_clr(T->row[0]+1, dim - 1);
	isl_seq_cpy(T->row[1], facet, dim);
	isl_seq_cpy(T->row[2], ridge, dim);
	T = isl_mat_right_inverse(T);
	set = isl_set_preimage(set, T);
	T = NULL;
	if (!set)
		goto error;
	lp = wrap_constraints(set);
	obj = isl_vec_alloc(ctx, 1 + dim*set->n);
	if (!obj)
		goto error;
	isl_int_set_si(obj->block.data[0], 0);
	for (i = 0; i < set->n; ++i) {
		isl_seq_clr(obj->block.data + 1 + dim*i, 2);
		isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
		isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
	}
	isl_int_init(num);
	isl_int_init(den);
	res = isl_basic_set_solve_lp(lp, 0,
			    obj->block.data, ctx->one, &num, &den, NULL);
	if (res == isl_lp_ok) {
		isl_int_neg(num, num);
		isl_seq_combine(facet, num, facet, den, ridge, dim);
		isl_seq_normalize(ctx, facet, dim);
	}
	isl_int_clear(num);
	isl_int_clear(den);
	isl_vec_free(obj);
	isl_basic_set_free(lp);
	isl_set_free(set);
	if (res == isl_lp_error)
		return NULL;
	isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded, 
		   return NULL);
	return facet;
error:
	isl_basic_set_free(lp);
	isl_mat_free(T);
	isl_set_free(set);
	return NULL;
}

/* Compute the constraint of a facet of "set".
 *
 * We first compute the intersection with a bounding constraint
 * that is orthogonal to one of the coordinate axes.
 * If the affine hull of this intersection has only one equality,
 * we have found a facet.
 * Otherwise, we wrap the current bounding constraint around
 * one of the equalities of the face (one that is not equal to
 * the current bounding constraint).
 * This process continues until we have found a facet.
 * The dimension of the intersection increases by at least
 * one on each iteration, so termination is guaranteed.
 */
static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
{
	struct isl_set *slice = NULL;
	struct isl_basic_set *face = NULL;
	int i;
	unsigned dim = isl_set_n_dim(set);
	int is_bound;
	isl_mat *bounds;

	isl_assert(set->ctx, set->n > 0, goto error);
	bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
	if (!bounds)
		return NULL;

	isl_seq_clr(bounds->row[0], dim);
	isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
	is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
	if (is_bound < 0)
		goto error;
	isl_assert(set->ctx, is_bound, goto error);
	isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
	bounds->n_row = 1;

	for (;;) {
		slice = isl_set_copy(set);
		slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
		face = isl_set_affine_hull(slice);
		if (!face)
			goto error;
		if (face->n_eq == 1) {
			isl_basic_set_free(face);
			break;
		}
		for (i = 0; i < face->n_eq; ++i)
			if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
			    !isl_seq_is_neg(bounds->row[0],
						face->eq[i], 1 + dim))
				break;
		isl_assert(set->ctx, i < face->n_eq, goto error);
		if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
			goto error;
		isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
		isl_basic_set_free(face);
	}

	return bounds;
error:
	isl_basic_set_free(face);
	isl_mat_free(bounds);
	return NULL;
}

/* Given the bounding constraint "c" of a facet of the convex hull of "set",
 * compute a hyperplane description of the facet, i.e., compute the facets
 * of the facet.
 *
 * We compute an affine transformation that transforms the constraint
 *
 *			  [ 1 ]
 *			c [ x ] = 0
 *
 * to the constraint
 *
 *			   z_1  = 0
 *
 * by computing the right inverse U of a matrix that starts with the rows
 *
 *			[ 1 0 ]
 *			[  c  ]
 *
 * Then
 *			[ 1 ]     [ 1 ]
 *			[ x ] = U [ z ]
 * and
 *			[ 1 ]     [ 1 ]
 *			[ z ] = Q [ x ]
 *
 * with Q = U^{-1}
 * Since z_1 is zero, we can drop this variable as well as the corresponding
 * column of U to obtain
 *
 *			[ 1 ]      [ 1  ]
 *			[ x ] = U' [ z' ]
 * and
 *			[ 1  ]      [ 1 ]
 *			[ z' ] = Q' [ x ]
 *
 * with Q' equal to Q, but without the corresponding row.
 * After computing the facets of the facet in the z' space,
 * we convert them back to the x space through Q.
 */
static struct isl_basic_set *compute_facet(struct isl_set *set, isl_int *c)
{
	struct isl_mat *m, *U, *Q;
	struct isl_basic_set *facet = NULL;
	struct isl_ctx *ctx;
	unsigned dim;

	ctx = set->ctx;
	set = isl_set_copy(set);
	dim = isl_set_n_dim(set);
	m = isl_mat_alloc(set->ctx, 2, 1 + dim);
	if (!m)
		goto error;
	isl_int_set_si(m->row[0][0], 1);
	isl_seq_clr(m->row[0]+1, dim);
	isl_seq_cpy(m->row[1], c, 1+dim);
	U = isl_mat_right_inverse(m);
	Q = isl_mat_right_inverse(isl_mat_copy(U));
	U = isl_mat_drop_cols(U, 1, 1);
	Q = isl_mat_drop_rows(Q, 1, 1);
	set = isl_set_preimage(set, U);
	facet = uset_convex_hull_wrap_bounded(set);
	facet = isl_basic_set_preimage(facet, Q);
	if (facet)
		isl_assert(ctx, facet->n_eq == 0, goto error);
	return facet;
error:
	isl_basic_set_free(facet);
	isl_set_free(set);
	return NULL;
}

/* Given an initial facet constraint, compute the remaining facets.
 * We do this by running through all facets found so far and computing
 * the adjacent facets through wrapping, adding those facets that we
 * hadn't already found before.
 *
 * For each facet we have found so far, we first compute its facets
 * in the resulting convex hull.  That is, we compute the ridges
 * of the resulting convex hull contained in the facet.
 * We also compute the corresponding facet in the current approximation
 * of the convex hull.  There is no need to wrap around the ridges
 * in this facet since that would result in a facet that is already
 * present in the current approximation.
 *
 * This function can still be significantly optimized by checking which of
 * the facets of the basic sets are also facets of the convex hull and
 * using all the facets so far to help in constructing the facets of the
 * facets
 * and/or
 * using the technique in section "3.1 Ridge Generation" of
 * "Extended Convex Hull" by Fukuda et al.
 */
static struct isl_basic_set *extend(struct isl_basic_set *hull,
	struct isl_set *set)
{
	int i, j, f;
	int k;
	struct isl_basic_set *facet = NULL;
	struct isl_basic_set *hull_facet = NULL;
	unsigned dim;

	if (!hull)
		return NULL;

	isl_assert(set->ctx, set->n > 0, goto error);

	dim = isl_set_n_dim(set);

	for (i = 0; i < hull->n_ineq; ++i) {
		facet = compute_facet(set, hull->ineq[i]);
		facet = isl_basic_set_add_equality(facet, hull->ineq[i]);
		facet = isl_basic_set_gauss(facet, NULL);
		facet = isl_basic_set_normalize_constraints(facet);
		hull_facet = isl_basic_set_copy(hull);
		hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]);
		hull_facet = isl_basic_set_gauss(hull_facet, NULL);
		hull_facet = isl_basic_set_normalize_constraints(hull_facet);
		if (!facet || !hull_facet)
			goto error;
		hull = isl_basic_set_cow(hull);
		hull = isl_basic_set_extend_dim(hull,
			isl_dim_copy(hull->dim), 0, 0, facet->n_ineq);
		if (!hull)
			goto error;
		for (j = 0; j < facet->n_ineq; ++j) {
			for (f = 0; f < hull_facet->n_ineq; ++f)
				if (isl_seq_eq(facet->ineq[j],
						hull_facet->ineq[f], 1 + dim))
					break;
			if (f < hull_facet->n_ineq)
				continue;
			k = isl_basic_set_alloc_inequality(hull);
			if (k < 0)
				goto error;
			isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
			if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
				goto error;
		}
		isl_basic_set_free(hull_facet);
		isl_basic_set_free(facet);
	}
	hull = isl_basic_set_simplify(hull);
	hull = isl_basic_set_finalize(hull);
	return hull;
error:
	isl_basic_set_free(hull_facet);
	isl_basic_set_free(facet);
	isl_basic_set_free(hull);
	return NULL;
}

/* Special case for computing the convex hull of a one dimensional set.
 * We simply collect the lower and upper bounds of each basic set
 * and the biggest of those.
 */
static struct isl_basic_set *convex_hull_1d(struct isl_set *set)
{
	struct isl_mat *c = NULL;
	isl_int *lower = NULL;
	isl_int *upper = NULL;
	int i, j, k;
	isl_int a, b;
	struct isl_basic_set *hull;

	for (i = 0; i < set->n; ++i) {
		set->p[i] = isl_basic_set_simplify(set->p[i]);
		if (!set->p[i])
			goto error;
	}
	set = isl_set_remove_empty_parts(set);
	if (!set)
		goto error;
	isl_assert(set->ctx, set->n > 0, goto error);
	c = isl_mat_alloc(set->ctx, 2, 2);
	if (!c)
		goto error;

	if (set->p[0]->n_eq > 0) {
		isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
		lower = c->row[0];
		upper = c->row[1];
		if (isl_int_is_pos(set->p[0]->eq[0][1])) {
			isl_seq_cpy(lower, set->p[0]->eq[0], 2);
			isl_seq_neg(upper, set->p[0]->eq[0], 2);
		} else {
			isl_seq_neg(lower, set->p[0]->eq[0], 2);
			isl_seq_cpy(upper, set->p[0]->eq[0], 2);
		}
	} else {
		for (j = 0; j < set->p[0]->n_ineq; ++j) {
			if (isl_int_is_pos(set->p[0]->ineq[j][1])) {
				lower = c->row[0];
				isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
			} else {
				upper = c->row[1];
				isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
			}
		}
	}

	isl_int_init(a);
	isl_int_init(b);
	for (i = 0; i < set->n; ++i) {
		struct isl_basic_set *bset = set->p[i];
		int has_lower = 0;
		int has_upper = 0;

		for (j = 0; j < bset->n_eq; ++j) {
			has_lower = 1;
			has_upper = 1;
			if (lower) {
				isl_int_mul(a, lower[0], bset->eq[j][1]);
				isl_int_mul(b, lower[1], bset->eq[j][0]);
				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
					isl_seq_cpy(lower, bset->eq[j], 2);
				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
					isl_seq_neg(lower, bset->eq[j], 2);
			}
			if (upper) {
				isl_int_mul(a, upper[0], bset->eq[j][1]);
				isl_int_mul(b, upper[1], bset->eq[j][0]);
				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))
					isl_seq_neg(upper, bset->eq[j], 2);
				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))
					isl_seq_cpy(upper, bset->eq[j], 2);
			}
		}
		for (j = 0; j < bset->n_ineq; ++j) {
			if (isl_int_is_pos(bset->ineq[j][1]))
				has_lower = 1;
			if (isl_int_is_neg(bset->ineq[j][1]))
				has_upper = 1;
			if (lower && isl_int_is_pos(bset->ineq[j][1])) {
				isl_int_mul(a, lower[0], bset->ineq[j][1]);
				isl_int_mul(b, lower[1], bset->ineq[j][0]);
				if (isl_int_lt(a, b))
					isl_seq_cpy(lower, bset->ineq[j], 2);
			}
			if (upper && isl_int_is_neg(bset->ineq[j][1])) {
				isl_int_mul(a, upper[0], bset->ineq[j][1]);
				isl_int_mul(b, upper[1], bset->ineq[j][0]);
				if (isl_int_gt(a, b))
					isl_seq_cpy(upper, bset->ineq[j], 2);
			}
		}
		if (!has_lower)
			lower = NULL;
		if (!has_upper)
			upper = NULL;
	}
	isl_int_clear(a);
	isl_int_clear(b);

	hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
	hull = isl_basic_set_set_rational(hull);
	if (!hull)
		goto error;
	if (lower) {
		k = isl_basic_set_alloc_inequality(hull);
		isl_seq_cpy(hull->ineq[k], lower, 2);
	}
	if (upper) {
		k = isl_basic_set_alloc_inequality(hull);
		isl_seq_cpy(hull->ineq[k], upper, 2);
	}
	hull = isl_basic_set_finalize(hull);
	isl_set_free(set);
	isl_mat_free(c);
	return hull;
error:
	isl_set_free(set);
	isl_mat_free(c);
	return NULL;
}

static struct isl_basic_set *convex_hull_0d(struct isl_set *set)
{
	struct isl_basic_set *convex_hull;

	if (!set)
		return NULL;

	if (isl_set_is_empty(set))
		convex_hull = isl_basic_set_empty(isl_dim_copy(set->dim));
	else
		convex_hull = isl_basic_set_universe(isl_dim_copy(set->dim));
	isl_set_free(set);
	return convex_hull;
}

/* Compute the convex hull of a pair of basic sets without any parameters or
 * integer divisions using Fourier-Motzkin elimination.
 * The convex hull is the set of all points that can be written as
 * the sum of points from both basic sets (in homogeneous coordinates).
 * We set up the constraints in a space with dimensions for each of
 * the three sets and then project out the dimensions corresponding
 * to the two original basic sets, retaining only those corresponding
 * to the convex hull.
 */
static struct isl_basic_set *convex_hull_pair_elim(struct isl_basic_set *bset1,
	struct isl_basic_set *bset2)
{
	int i, j, k;
{'fraction_non_alphanumeric': 0.1409131115193577, 'fraction_numerical': 0.013109251321131297, 'mean_word_length': 3.899655517279698, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 423, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2593841', 'n_tokens_mistral': 20103, 'n_tokens_neox': 18520, 'n_words': 4716}