/*

 * Copyright 2008-2009 Katholieke Universiteit Leuven

 * Copyright 2014      INRIA Rocquencourt

 *

 * Use of this software is governed by the MIT license

 *

 * Written by Sven Verdoolaege, K.U.Leuven, Departement

 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium

 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,

 * B.P. 105 - 78153 Le Chesnay, France

 */

#include <isl_ctx_private.h>

#include <isl_map_private.h>

#include <isl_lp_private.h>

#include <isl/map.h>

#include <isl_mat_private.h>

#include <isl_vec_private.h>

#include <isl/set.h>

#include <isl_seq.h>

#include <isl_options_private.h>

#include "isl_equalities.h"

#include "isl_tab.h"

#include <isl_sort.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, 0);

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

}

/* Remove redundant constraints in each of the basic maps.

 */

__isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)

{

	return isl_map_inline_foreach_basic_map(map,

					    &isl_basic_map_remove_redundancies);

}

__isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)

{

	return isl_map_remove_redundancies(set);

}

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

}

__isl_give isl_map *isl_map_set_rational(__isl_take isl_map *map)

{

	int i;

	map = isl_map_cow(map);

	if (!map)

		return NULL;

	for (i = 0; i < map->n; ++i) {

		map->p[i] = isl_basic_map_set_rational(map->p[i]);

		if (!map->p[i])

			goto error;

	}

	return map;

error:

	isl_map_free(map);

	return NULL;

}

__isl_give isl_set *isl_set_set_rational(__isl_take isl_set *set)

{

	return isl_map_set_rational(set);

}

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 = NULL;

	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_space(hull,

			isl_space_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_space_copy(set->dim));

	else

		convex_hull = isl_basic_set_universe(isl_space_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;

	struct isl_basic_set *bset[2];

	struct isl_basic_set *hull = NULL;

	unsigned dim;

	if (!bset1 || !bset2)

		goto error;

	dim = isl_basic_set_n_dim(bset1);

	hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,

				1 + dim + bset1->n_eq + bset2->n_eq,

				2 + bset1->n_ineq + bset2->n_ineq);

	bset[0] = bset1;

	bset[1] = bset2;

	for (i = 0; i < 2; ++i) {

		for (j = 0; j < bset[i]->n_eq; ++j) {

			k = isl_basic_set_alloc_equality(hull);

			if (k < 0)

				goto error;

			isl_seq_clr(hull->eq[k], (i+1) * (1+dim));

			isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));

			isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],

					1+dim);

		}

		for (j = 0; j < bset[i]->n_ineq; ++j) {

			k = isl_basic_set_alloc_inequality(hull);

			if (k < 0)

				goto error;

			isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));

			isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));

			isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),

					bset[i]->ineq[j], 1+dim);

		}

		k = isl_basic_set_alloc_inequality(hull);

		if (k < 0)

			goto error;

		isl_seq_clr(hull->ineq[k], 1+2+3*dim);

		isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);

	}

	for (j = 0; j < 1+dim; ++j) {

		k = isl_basic_set_alloc_equality(hull);

		if (k < 0)

			goto error;