brintos / llvm-project-archived public Read only

1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 * Copyright 2014      INRIA Rocquencourt4 *5 * Use of this software is governed by the MIT license6 *7 * Written by Sven Verdoolaege, K.U.Leuven, Departement8 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium9 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,10 * B.P. 105 - 78153 Le Chesnay, France11 */12 13#include <isl_ctx_private.h>14#include <isl_map_private.h>15#include <isl_lp_private.h>16#include <isl/map.h>17#include <isl_mat_private.h>18#include <isl_vec_private.h>19#include <isl/set.h>20#include <isl_seq.h>21#include <isl_options_private.h>22#include "isl_equalities.h"23#include "isl_tab.h"24#include <isl_sort.h>25 26#include <bset_to_bmap.c>27#include <bset_from_bmap.c>28#include <set_to_map.c>29 30static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(31	__isl_take isl_set *set);32 33/* Remove redundant34 * constraints.  If the minimal value along the normal of a constraint35 * is the same if the constraint is removed, then the constraint is redundant.36 *37 * Since some constraints may be mutually redundant, sort the constraints38 * first such that constraints that involve existentially quantified39 * variables are considered for removal before those that do not.40 * The sorting is also needed for the use in map_simple_hull.41 *42 * Note that isl_tab_detect_implicit_equalities may also end up43 * marking some constraints as redundant.  Make sure the constraints44 * are preserved and undo those marking such that isl_tab_detect_redundant45 * can consider the constraints in the sorted order.46 *47 * Alternatively, we could have intersected the basic map with the48 * corresponding equality and then checked if the dimension was that49 * of a facet.50 */51__isl_give isl_basic_map *isl_basic_map_remove_redundancies(52	__isl_take isl_basic_map *bmap)53{54	struct isl_tab *tab;55 56	if (!bmap)57		return NULL;58 59	bmap = isl_basic_map_gauss(bmap, NULL);60	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))61		return bmap;62	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT))63		return bmap;64	if (bmap->n_ineq <= 1)65		return bmap;66 67	bmap = isl_basic_map_sort_constraints(bmap);68	tab = isl_tab_from_basic_map(bmap, 0);69	if (!tab)70		goto error;71	tab->preserve = 1;72	if (isl_tab_detect_implicit_equalities(tab) < 0)73		goto error;74	if (isl_tab_restore_redundant(tab) < 0)75		goto error;76	tab->preserve = 0;77	if (isl_tab_detect_redundant(tab) < 0)78		goto error;79	bmap = isl_basic_map_update_from_tab(bmap, tab);80	isl_tab_free(tab);81	if (!bmap)82		return NULL;83	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);84	ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);85	return bmap;86error:87	isl_tab_free(tab);88	isl_basic_map_free(bmap);89	return NULL;90}91 92__isl_give isl_basic_set *isl_basic_set_remove_redundancies(93	__isl_take isl_basic_set *bset)94{95	return bset_from_bmap(96		isl_basic_map_remove_redundancies(bset_to_bmap(bset)));97}98 99/* Remove redundant constraints in each of the basic maps.100 */101__isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)102{103	return isl_map_inline_foreach_basic_map(map,104					    &isl_basic_map_remove_redundancies);105}106 107__isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)108{109	return isl_map_remove_redundancies(set);110}111 112/* Check if the set set is bound in the direction of the affine113 * constraint c and if so, set the constant term such that the114 * resulting constraint is a bounding constraint for the set.115 */116static isl_bool uset_is_bound(__isl_keep isl_set *set, isl_int *c, unsigned len)117{118	int first;119	int j;120	isl_int opt;121	isl_int opt_denom;122 123	isl_int_init(opt);124	isl_int_init(opt_denom);125	first = 1;126	for (j = 0; j < set->n; ++j) {127		enum isl_lp_result res;128 129		if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))130			continue;131 132		res = isl_basic_set_solve_lp(set->p[j],133				0, c, set->ctx->one, &opt, &opt_denom, NULL);134		if (res == isl_lp_unbounded)135			break;136		if (res == isl_lp_error)137			goto error;138		if (res == isl_lp_empty) {139			set->p[j] = isl_basic_set_set_to_empty(set->p[j]);140			if (!set->p[j])141				goto error;142			continue;143		}144		if (first || isl_int_is_neg(opt)) {145			if (!isl_int_is_one(opt_denom))146				isl_seq_scale(c, c, opt_denom, len);147			isl_int_sub(c[0], c[0], opt);148		}149		first = 0;150	}151	isl_int_clear(opt);152	isl_int_clear(opt_denom);153	return isl_bool_ok(j >= set->n);154error:155	isl_int_clear(opt);156	isl_int_clear(opt_denom);157	return isl_bool_error;158}159 160static __isl_give isl_set *isl_set_add_basic_set_equality(161	__isl_take isl_set *set, isl_int *c)162{163	int i;164 165	set = isl_set_cow(set);166	if (!set)167		return NULL;168	for (i = 0; i < set->n; ++i) {169		set->p[i] = isl_basic_set_add_eq(set->p[i], c);170		if (!set->p[i])171			goto error;172	}173	return set;174error:175	isl_set_free(set);176	return NULL;177}178 179/* Given a union of basic sets, construct the constraints for wrapping180 * a facet around one of its ridges.181 * In particular, if each of n the d-dimensional basic sets i in "set"182 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0183 * and is defined by the constraints184 *				    [ 1 ]185 *				A_i [ x ]  >= 0186 *187 * then the resulting set is of dimension n*(1+d) and has as constraints188 *189 *				    [ a_i ]190 *				A_i [ x_i ] >= 0191 *192 *				      a_i   >= 0193 *194 *			\sum_i x_{i,1} = 1195 */196static __isl_give isl_basic_set *wrap_constraints(__isl_keep isl_set *set)197{198	struct isl_basic_set *lp;199	unsigned n_eq;200	unsigned n_ineq;201	int i, j, k;202	isl_size dim, lp_dim;203 204	dim = isl_set_dim(set, isl_dim_set);205	if (dim < 0)206		return NULL;207 208	dim += 1;209	n_eq = 1;210	n_ineq = set->n;211	for (i = 0; i < set->n; ++i) {212		n_eq += set->p[i]->n_eq;213		n_ineq += set->p[i]->n_ineq;214	}215	lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);216	lp = isl_basic_set_set_rational(lp);217	if (!lp)218		return NULL;219	lp_dim = isl_basic_set_dim(lp, isl_dim_set);220	if (lp_dim < 0)221		return isl_basic_set_free(lp);222	k = isl_basic_set_alloc_equality(lp);223	isl_int_set_si(lp->eq[k][0], -1);224	for (i = 0; i < set->n; ++i) {225		isl_int_set_si(lp->eq[k][1+dim*i], 0);226		isl_int_set_si(lp->eq[k][1+dim*i+1], 1);227		isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);228	}229	for (i = 0; i < set->n; ++i) {230		k = isl_basic_set_alloc_inequality(lp);231		isl_seq_clr(lp->ineq[k], 1+lp_dim);232		isl_int_set_si(lp->ineq[k][1+dim*i], 1);233 234		for (j = 0; j < set->p[i]->n_eq; ++j) {235			k = isl_basic_set_alloc_equality(lp);236			isl_seq_clr(lp->eq[k], 1+dim*i);237			isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);238			isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));239		}240 241		for (j = 0; j < set->p[i]->n_ineq; ++j) {242			k = isl_basic_set_alloc_inequality(lp);243			isl_seq_clr(lp->ineq[k], 1+dim*i);244			isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);245			isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));246		}247	}248	return lp;249}250 251/* Given a facet "facet" of the convex hull of "set" and a facet "ridge"252 * of that facet, compute the other facet of the convex hull that contains253 * the ridge.254 *255 * We first transform the set such that the facet constraint becomes256 *257 *			x_1 >= 0258 *259 * I.e., the facet lies in260 *261 *			x_1 = 0262 *263 * and on that facet, the constraint that defines the ridge is264 *265 *			x_2 >= 0266 *267 * (This transformation is not strictly needed, all that is needed is268 * that the ridge contains the origin.)269 *270 * Since the ridge contains the origin, the cone of the convex hull271 * will be of the form272 *273 *			x_1 >= 0274 *			x_2 >= a x_1275 *276 * with this second constraint defining the new facet.277 * The constant a is obtained by settting x_1 in the cone of the278 * convex hull to 1 and minimizing x_2.279 * Now, each element in the cone of the convex hull is the sum280 * of elements in the cones of the basic sets.281 * If a_i is the dilation factor of basic set i, then the problem282 * we need to solve is283 *284 *			min \sum_i x_{i,2}285 *			st286 *				\sum_i x_{i,1} = 1287 *				    a_i   >= 0288 *				  [ a_i ]289 *				A [ x_i ] >= 0290 *291 * with292 *				    [  1  ]293 *				A_i [ x_i ] >= 0294 *295 * the constraints of each (transformed) basic set.296 * If a = n/d, then the constraint defining the new facet (in the transformed297 * space) is298 *299 *			-n x_1 + d x_2 >= 0300 *301 * In the original space, we need to take the same combination of the302 * corresponding constraints "facet" and "ridge".303 *304 * If a = -infty = "-1/0", then we just return the original facet constraint.305 * This means that the facet is unbounded, but has a bounded intersection306 * with the union of sets.307 */308isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,309	isl_int *facet, isl_int *ridge)310{311	int i;312	isl_ctx *ctx;313	struct isl_mat *T = NULL;314	struct isl_basic_set *lp = NULL;315	struct isl_vec *obj;316	enum isl_lp_result res;317	isl_int num, den;318	isl_size dim;319 320	dim = isl_set_dim(set, isl_dim_set);321	if (dim < 0)322		return NULL;323	ctx = set->ctx;324	set = isl_set_copy(set);325	set = isl_set_set_rational(set);326 327	dim += 1;328	T = isl_mat_alloc(ctx, 3, dim);329	if (!T)330		goto error;331	isl_int_set_si(T->row[0][0], 1);332	isl_seq_clr(T->row[0]+1, dim - 1);333	isl_seq_cpy(T->row[1], facet, dim);334	isl_seq_cpy(T->row[2], ridge, dim);335	T = isl_mat_right_inverse(T);336	set = isl_set_preimage(set, T);337	T = NULL;338	if (!set)339		goto error;340	lp = wrap_constraints(set);341	obj = isl_vec_alloc(ctx, 1 + dim*set->n);342	if (!obj)343		goto error;344	isl_int_set_si(obj->block.data[0], 0);345	for (i = 0; i < set->n; ++i) {346		isl_seq_clr(obj->block.data + 1 + dim*i, 2);347		isl_int_set_si(obj->block.data[1 + dim*i+2], 1);348		isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);349	}350	isl_int_init(num);351	isl_int_init(den);352	res = isl_basic_set_solve_lp(lp, 0,353			    obj->block.data, ctx->one, &num, &den, NULL);354	if (res == isl_lp_ok) {355		isl_int_neg(num, num);356		isl_seq_combine(facet, num, facet, den, ridge, dim);357		isl_seq_normalize(ctx, facet, dim);358	}359	isl_int_clear(num);360	isl_int_clear(den);361	isl_vec_free(obj);362	isl_basic_set_free(lp);363	isl_set_free(set);364	if (res == isl_lp_error)365		return NULL;366	isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded, 367		   return NULL);368	return facet;369error:370	isl_basic_set_free(lp);371	isl_mat_free(T);372	isl_set_free(set);373	return NULL;374}375 376/* Compute the constraint of a facet of "set".377 *378 * We first compute the intersection with a bounding constraint379 * that is orthogonal to one of the coordinate axes.380 * If the affine hull of this intersection has only one equality,381 * we have found a facet.382 * Otherwise, we wrap the current bounding constraint around383 * one of the equalities of the face (one that is not equal to384 * the current bounding constraint).385 * This process continues until we have found a facet.386 * The dimension of the intersection increases by at least387 * one on each iteration, so termination is guaranteed.388 */389static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)390{391	struct isl_set *slice = NULL;392	struct isl_basic_set *face = NULL;393	int i;394	isl_size dim = isl_set_dim(set, isl_dim_set);395	isl_bool is_bound;396	isl_mat *bounds = NULL;397 398	if (dim < 0)399		return NULL;400	isl_assert(set->ctx, set->n > 0, goto error);401	bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);402	if (!bounds)403		return NULL;404 405	isl_seq_clr(bounds->row[0], dim);406	isl_int_set_si(bounds->row[0][1 + dim - 1], 1);407	is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);408	if (is_bound < 0)409		goto error;410	isl_assert(set->ctx, is_bound, goto error);411	isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);412	bounds->n_row = 1;413 414	for (;;) {415		slice = isl_set_copy(set);416		slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);417		face = isl_set_affine_hull(slice);418		if (!face)419			goto error;420		if (face->n_eq == 1) {421			isl_basic_set_free(face);422			break;423		}424		for (i = 0; i < face->n_eq; ++i)425			if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&426			    !isl_seq_is_neg(bounds->row[0],427						face->eq[i], 1 + dim))428				break;429		isl_assert(set->ctx, i < face->n_eq, goto error);430		if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))431			goto error;432		isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);433		isl_basic_set_free(face);434	}435 436	return bounds;437error:438	isl_basic_set_free(face);439	isl_mat_free(bounds);440	return NULL;441}442 443/* Given the bounding constraint "c" of a facet of the convex hull of "set",444 * compute a hyperplane description of the facet, i.e., compute the facets445 * of the facet.446 *447 * We compute an affine transformation that transforms the constraint448 *449 *			  [ 1 ]450 *			c [ x ] = 0451 *452 * to the constraint453 *454 *			   z_1  = 0455 *456 * by computing the right inverse U of a matrix that starts with the rows457 *458 *			[ 1 0 ]459 *			[  c  ]460 *461 * Then462 *			[ 1 ]     [ 1 ]463 *			[ x ] = U [ z ]464 * and465 *			[ 1 ]     [ 1 ]466 *			[ z ] = Q [ x ]467 *468 * with Q = U^{-1}469 * Since z_1 is zero, we can drop this variable as well as the corresponding470 * column of U to obtain471 *472 *			[ 1 ]      [ 1  ]473 *			[ x ] = U' [ z' ]474 * and475 *			[ 1  ]      [ 1 ]476 *			[ z' ] = Q' [ x ]477 *478 * with Q' equal to Q, but without the corresponding row.479 * After computing the facets of the facet in the z' space,480 * we convert them back to the x space through Q.481 */482static __isl_give isl_basic_set *compute_facet(__isl_keep isl_set *set,483	isl_int *c)484{485	struct isl_mat *m, *U, *Q;486	struct isl_basic_set *facet = NULL;487	struct isl_ctx *ctx;488	isl_size dim;489 490	dim = isl_set_dim(set, isl_dim_set);491	if (dim < 0)492		return NULL;493	ctx = set->ctx;494	set = isl_set_copy(set);495	m = isl_mat_alloc(set->ctx, 2, 1 + dim);496	if (!m)497		goto error;498	isl_int_set_si(m->row[0][0], 1);499	isl_seq_clr(m->row[0]+1, dim);500	isl_seq_cpy(m->row[1], c, 1+dim);501	U = isl_mat_right_inverse(m);502	Q = isl_mat_right_inverse(isl_mat_copy(U));503	U = isl_mat_drop_cols(U, 1, 1);504	Q = isl_mat_drop_rows(Q, 1, 1);505	set = isl_set_preimage(set, U);506	facet = uset_convex_hull_wrap_bounded(set);507	facet = isl_basic_set_preimage(facet, Q);508	if (facet && facet->n_eq != 0)509		isl_die(ctx, isl_error_internal, "unexpected equality",510			return isl_basic_set_free(facet));511	return facet;512error:513	isl_basic_set_free(facet);514	isl_set_free(set);515	return NULL;516}517 518/* Given an initial facet constraint, compute the remaining facets.519 * We do this by running through all facets found so far and computing520 * the adjacent facets through wrapping, adding those facets that we521 * hadn't already found before.522 *523 * For each facet we have found so far, we first compute its facets524 * in the resulting convex hull.  That is, we compute the ridges525 * of the resulting convex hull contained in the facet.526 * We also compute the corresponding facet in the current approximation527 * of the convex hull.  There is no need to wrap around the ridges528 * in this facet since that would result in a facet that is already529 * present in the current approximation.530 *531 * This function can still be significantly optimized by checking which of532 * the facets of the basic sets are also facets of the convex hull and533 * using all the facets so far to help in constructing the facets of the534 * facets535 * and/or536 * using the technique in section "3.1 Ridge Generation" of537 * "Extended Convex Hull" by Fukuda et al.538 */539static __isl_give isl_basic_set *extend(__isl_take isl_basic_set *hull,540	__isl_keep isl_set *set)541{542	int i, j, f;543	int k;544	struct isl_basic_set *facet = NULL;545	struct isl_basic_set *hull_facet = NULL;546	isl_size dim;547 548	dim = isl_set_dim(set, isl_dim_set);549	if (dim < 0 || !hull)550		return isl_basic_set_free(hull);551 552	isl_assert(set->ctx, set->n > 0, goto error);553 554	for (i = 0; i < hull->n_ineq; ++i) {555		facet = compute_facet(set, hull->ineq[i]);556		facet = isl_basic_set_add_eq(facet, hull->ineq[i]);557		facet = isl_basic_set_gauss(facet, NULL);558		facet = isl_basic_set_normalize_constraints(facet);559		hull_facet = isl_basic_set_copy(hull);560		hull_facet = isl_basic_set_add_eq(hull_facet, hull->ineq[i]);561		hull_facet = isl_basic_set_gauss(hull_facet, NULL);562		hull_facet = isl_basic_set_normalize_constraints(hull_facet);563		if (!facet || !hull_facet)564			goto error;565		hull = isl_basic_set_cow(hull);566		hull = isl_basic_set_extend(hull, 0, 0, facet->n_ineq);567		if (!hull)568			goto error;569		for (j = 0; j < facet->n_ineq; ++j) {570			for (f = 0; f < hull_facet->n_ineq; ++f)571				if (isl_seq_eq(facet->ineq[j],572						hull_facet->ineq[f], 1 + dim))573					break;574			if (f < hull_facet->n_ineq)575				continue;576			k = isl_basic_set_alloc_inequality(hull);577			if (k < 0)578				goto error;579			isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);580			if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))581				goto error;582		}583		isl_basic_set_free(hull_facet);584		isl_basic_set_free(facet);585	}586	hull = isl_basic_set_simplify(hull);587	hull = isl_basic_set_finalize(hull);588	return hull;589error:590	isl_basic_set_free(hull_facet);591	isl_basic_set_free(facet);592	isl_basic_set_free(hull);593	return NULL;594}595 596/* Special case for computing the convex hull of a one dimensional set.597 * We simply collect the lower and upper bounds of each basic set598 * and the biggest of those.599 */600static __isl_give isl_basic_set *convex_hull_1d(__isl_take isl_set *set)601{602	struct isl_mat *c = NULL;603	isl_int *lower = NULL;604	isl_int *upper = NULL;605	int i, j, k;606	isl_int a, b;607	struct isl_basic_set *hull;608 609	for (i = 0; i < set->n; ++i) {610		set->p[i] = isl_basic_set_simplify(set->p[i]);611		if (!set->p[i])612			goto error;613	}614	set = isl_set_remove_empty_parts(set);615	if (!set)616		goto error;617	isl_assert(set->ctx, set->n > 0, goto error);618	c = isl_mat_alloc(set->ctx, 2, 2);619	if (!c)620		goto error;621 622	if (set->p[0]->n_eq > 0) {623		isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);624		lower = c->row[0];625		upper = c->row[1];626		if (isl_int_is_pos(set->p[0]->eq[0][1])) {627			isl_seq_cpy(lower, set->p[0]->eq[0], 2);628			isl_seq_neg(upper, set->p[0]->eq[0], 2);629		} else {630			isl_seq_neg(lower, set->p[0]->eq[0], 2);631			isl_seq_cpy(upper, set->p[0]->eq[0], 2);632		}633	} else {634		for (j = 0; j < set->p[0]->n_ineq; ++j) {635			if (isl_int_is_pos(set->p[0]->ineq[j][1])) {636				lower = c->row[0];637				isl_seq_cpy(lower, set->p[0]->ineq[j], 2);638			} else {639				upper = c->row[1];640				isl_seq_cpy(upper, set->p[0]->ineq[j], 2);641			}642		}643	}644 645	isl_int_init(a);646	isl_int_init(b);647	for (i = 0; i < set->n; ++i) {648		struct isl_basic_set *bset = set->p[i];649		int has_lower = 0;650		int has_upper = 0;651 652		for (j = 0; j < bset->n_eq; ++j) {653			has_lower = 1;654			has_upper = 1;655			if (lower) {656				isl_int_mul(a, lower[0], bset->eq[j][1]);657				isl_int_mul(b, lower[1], bset->eq[j][0]);658				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))659					isl_seq_cpy(lower, bset->eq[j], 2);660				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))661					isl_seq_neg(lower, bset->eq[j], 2);662			}663			if (upper) {664				isl_int_mul(a, upper[0], bset->eq[j][1]);665				isl_int_mul(b, upper[1], bset->eq[j][0]);666				if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1]))667					isl_seq_neg(upper, bset->eq[j], 2);668				if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1]))669					isl_seq_cpy(upper, bset->eq[j], 2);670			}671		}672		for (j = 0; j < bset->n_ineq; ++j) {673			if (isl_int_is_pos(bset->ineq[j][1]))674				has_lower = 1;675			if (isl_int_is_neg(bset->ineq[j][1]))676				has_upper = 1;677			if (lower && isl_int_is_pos(bset->ineq[j][1])) {678				isl_int_mul(a, lower[0], bset->ineq[j][1]);679				isl_int_mul(b, lower[1], bset->ineq[j][0]);680				if (isl_int_lt(a, b))681					isl_seq_cpy(lower, bset->ineq[j], 2);682			}683			if (upper && isl_int_is_neg(bset->ineq[j][1])) {684				isl_int_mul(a, upper[0], bset->ineq[j][1]);685				isl_int_mul(b, upper[1], bset->ineq[j][0]);686				if (isl_int_gt(a, b))687					isl_seq_cpy(upper, bset->ineq[j], 2);688			}689		}690		if (!has_lower)691			lower = NULL;692		if (!has_upper)693			upper = NULL;694	}695	isl_int_clear(a);696	isl_int_clear(b);697 698	hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);699	hull = isl_basic_set_set_rational(hull);700	if (!hull)701		goto error;702	if (lower) {703		k = isl_basic_set_alloc_inequality(hull);704		isl_seq_cpy(hull->ineq[k], lower, 2);705	}706	if (upper) {707		k = isl_basic_set_alloc_inequality(hull);708		isl_seq_cpy(hull->ineq[k], upper, 2);709	}710	hull = isl_basic_set_finalize(hull);711	isl_set_free(set);712	isl_mat_free(c);713	return hull;714error:715	isl_set_free(set);716	isl_mat_free(c);717	return NULL;718}719 720static __isl_give isl_basic_set *convex_hull_0d(__isl_take isl_set *set)721{722	struct isl_basic_set *convex_hull;723 724	if (!set)725		return NULL;726 727	if (isl_set_is_empty(set))728		convex_hull = isl_basic_set_empty(isl_space_copy(set->dim));729	else730		convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));731	isl_set_free(set);732	return convex_hull;733}734 735/* Compute the convex hull of a pair of basic sets without any parameters or736 * integer divisions using Fourier-Motzkin elimination.737 * The convex hull is the set of all points that can be written as738 * the sum of points from both basic sets (in homogeneous coordinates).739 * We set up the constraints in a space with dimensions for each of740 * the three sets and then project out the dimensions corresponding741 * to the two original basic sets, retaining only those corresponding742 * to the convex hull.743 */744static __isl_give isl_basic_set *convex_hull_pair_elim(745	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)746{747	int i, j, k;748	struct isl_basic_set *bset[2];749	struct isl_basic_set *hull = NULL;750	isl_size dim;751 752	dim = isl_basic_set_dim(bset1, isl_dim_set);753	if (dim < 0 || !bset2)754		goto error;755 756	hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,757				1 + dim + bset1->n_eq + bset2->n_eq,758				2 + bset1->n_ineq + bset2->n_ineq);759	bset[0] = bset1;760	bset[1] = bset2;761	for (i = 0; i < 2; ++i) {762		for (j = 0; j < bset[i]->n_eq; ++j) {763			k = isl_basic_set_alloc_equality(hull);764			if (k < 0)765				goto error;766			isl_seq_clr(hull->eq[k], (i+1) * (1+dim));767			isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));768			isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],769					1+dim);770		}771		for (j = 0; j < bset[i]->n_ineq; ++j) {772			k = isl_basic_set_alloc_inequality(hull);773			if (k < 0)774				goto error;775			isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));776			isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));777			isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),778					bset[i]->ineq[j], 1+dim);779		}780		k = isl_basic_set_alloc_inequality(hull);781		if (k < 0)782			goto error;783		isl_seq_clr(hull->ineq[k], 1+2+3*dim);784		isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);785	}786	for (j = 0; j < 1+dim; ++j) {787		k = isl_basic_set_alloc_equality(hull);788		if (k < 0)789			goto error;790		isl_seq_clr(hull->eq[k], 1+2+3*dim);791		isl_int_set_si(hull->eq[k][j], -1);792		isl_int_set_si(hull->eq[k][1+dim+j], 1);793		isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);794	}795	hull = isl_basic_set_set_rational(hull);796	hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim));797	hull = isl_basic_set_remove_redundancies(hull);798	isl_basic_set_free(bset1);799	isl_basic_set_free(bset2);800	return hull;801error:802	isl_basic_set_free(bset1);803	isl_basic_set_free(bset2);804	isl_basic_set_free(hull);805	return NULL;806}807 808/* Is the set bounded for each value of the parameters?809 */810isl_bool isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)811{812	struct isl_tab *tab;813	isl_bool bounded;814 815	if (!bset)816		return isl_bool_error;817	if (isl_basic_set_plain_is_empty(bset))818		return isl_bool_true;819 820	tab = isl_tab_from_recession_cone(bset, 1);821	bounded = isl_tab_cone_is_bounded(tab);822	isl_tab_free(tab);823	return bounded;824}825 826/* Is the image bounded for each value of the parameters and827 * the domain variables?828 */829isl_bool isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)830{831	isl_size nparam = isl_basic_map_dim(bmap, isl_dim_param);832	isl_size n_in = isl_basic_map_dim(bmap, isl_dim_in);833	isl_bool bounded;834 835	if (nparam < 0 || n_in < 0)836		return isl_bool_error;837 838	bmap = isl_basic_map_copy(bmap);839	bmap = isl_basic_map_cow(bmap);840	bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,841					isl_dim_in, 0, n_in);842	bounded = isl_basic_set_is_bounded(bset_from_bmap(bmap));843	isl_basic_map_free(bmap);844 845	return bounded;846}847 848/* Is the set bounded for each value of the parameters?849 */850isl_bool isl_set_is_bounded(__isl_keep isl_set *set)851{852	int i;853 854	if (!set)855		return isl_bool_error;856 857	for (i = 0; i < set->n; ++i) {858		isl_bool bounded = isl_basic_set_is_bounded(set->p[i]);859		if (!bounded || bounded < 0)860			return bounded;861	}862	return isl_bool_true;863}864 865/* Compute the lineality space of the convex hull of bset1 and bset2.866 *867 * We first compute the intersection of the recession cone of bset1868 * with the negative of the recession cone of bset2 and then compute869 * the linear hull of the resulting cone.870 */871static __isl_give isl_basic_set *induced_lineality_space(872	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)873{874	int i, k;875	struct isl_basic_set *lin = NULL;876	isl_size dim;877 878	dim = isl_basic_set_dim(bset1, isl_dim_all);879	if (dim < 0 || !bset2)880		goto error;881 882	lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0,883					bset1->n_eq + bset2->n_eq,884					bset1->n_ineq + bset2->n_ineq);885	lin = isl_basic_set_set_rational(lin);886	if (!lin)887		goto error;888	for (i = 0; i < bset1->n_eq; ++i) {889		k = isl_basic_set_alloc_equality(lin);890		if (k < 0)891			goto error;892		isl_int_set_si(lin->eq[k][0], 0);893		isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim);894	}895	for (i = 0; i < bset1->n_ineq; ++i) {896		k = isl_basic_set_alloc_inequality(lin);897		if (k < 0)898			goto error;899		isl_int_set_si(lin->ineq[k][0], 0);900		isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim);901	}902	for (i = 0; i < bset2->n_eq; ++i) {903		k = isl_basic_set_alloc_equality(lin);904		if (k < 0)905			goto error;906		isl_int_set_si(lin->eq[k][0], 0);907		isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim);908	}909	for (i = 0; i < bset2->n_ineq; ++i) {910		k = isl_basic_set_alloc_inequality(lin);911		if (k < 0)912			goto error;913		isl_int_set_si(lin->ineq[k][0], 0);914		isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim);915	}916 917	isl_basic_set_free(bset1);918	isl_basic_set_free(bset2);919	return isl_basic_set_affine_hull(lin);920error:921	isl_basic_set_free(lin);922	isl_basic_set_free(bset1);923	isl_basic_set_free(bset2);924	return NULL;925}926 927static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set);928 929/* Given a set and a linear space "lin" of dimension n > 0,930 * project the linear space from the set, compute the convex hull931 * and then map the set back to the original space.932 *933 * Let934 *935 *	M x = 0936 *937 * describe the linear space.  We first compute the Hermite normal938 * form H = M U of M = H Q, to obtain939 *940 *	H Q x = 0941 *942 * The last n rows of H will be zero, so the last n variables of x' = Q x943 * are the one we want to project out.  We do this by transforming each944 * basic set A x >= b to A U x' >= b and then removing the last n dimensions.945 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b',946 * we transform the hull back to the original space as A' Q_1 x >= b',947 * with Q_1 all but the last n rows of Q.948 */949static __isl_give isl_basic_set *modulo_lineality(__isl_take isl_set *set,950	__isl_take isl_basic_set *lin)951{952	isl_size total = isl_basic_set_dim(lin, isl_dim_all);953	unsigned lin_dim;954	struct isl_basic_set *hull;955	struct isl_mat *M, *U, *Q;956 957	if (!set || total < 0)958		goto error;959	lin_dim = total - lin->n_eq;960	M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total);961	M = isl_mat_left_hermite(M, 0, &U, &Q);962	if (!M)963		goto error;964	isl_mat_free(M);965	isl_basic_set_free(lin);966 967	Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim);968 969	U = isl_mat_lin_to_aff(U);970	Q = isl_mat_lin_to_aff(Q);971 972	set = isl_set_preimage(set, U);973	set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim);974	hull = uset_convex_hull(set);975	hull = isl_basic_set_preimage(hull, Q);976 977	return hull;978error:979	isl_basic_set_free(lin);980	isl_set_free(set);981	return NULL;982}983 984/* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,985 * set up an LP for solving986 *987 *	\sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}988 *989 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0990 * The next \alpha{ij} correspond to the equalities and come in pairs.991 * The final \alpha{ij} correspond to the inequalities.992 */993static __isl_give isl_basic_set *valid_direction_lp(994	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)995{996	isl_space *space;997	struct isl_basic_set *lp;998	unsigned d;999	int n;1000	int i, j, k;1001	isl_size total;1002 1003	total = isl_basic_set_dim(bset1, isl_dim_all);1004	if (total < 0 || !bset2)1005		goto error;1006	d = 1 + total;1007	n = 2 +1008	    2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;1009	space = isl_space_set_alloc(bset1->ctx, 0, n);1010	lp = isl_basic_set_alloc_space(space, 0, d, n);1011	if (!lp)1012		goto error;1013	for (i = 0; i < n; ++i) {1014		k = isl_basic_set_alloc_inequality(lp);1015		if (k < 0)1016			goto error;1017		isl_seq_clr(lp->ineq[k] + 1, n);1018		isl_int_set_si(lp->ineq[k][0], -1);1019		isl_int_set_si(lp->ineq[k][1 + i], 1);1020	}1021	for (i = 0; i < d; ++i) {1022		k = isl_basic_set_alloc_equality(lp);1023		if (k < 0)1024			goto error;1025		n = 0;1026		isl_int_set_si(lp->eq[k][n], 0); n++;1027		/* positivity constraint 1 >= 0 */1028		isl_int_set_si(lp->eq[k][n], i == 0); n++;1029		for (j = 0; j < bset1->n_eq; ++j) {1030			isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;1031			isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;1032		}1033		for (j = 0; j < bset1->n_ineq; ++j) {1034			isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;1035		}1036		/* positivity constraint 1 >= 0 */1037		isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;1038		for (j = 0; j < bset2->n_eq; ++j) {1039			isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;1040			isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;1041		}1042		for (j = 0; j < bset2->n_ineq; ++j) {1043			isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;1044		}1045	}1046	lp = isl_basic_set_gauss(lp, NULL);1047	isl_basic_set_free(bset1);1048	isl_basic_set_free(bset2);1049	return lp;1050error:1051	isl_basic_set_free(bset1);1052	isl_basic_set_free(bset2);1053	return NULL;1054}1055 1056/* Compute a vector s in the homogeneous space such that <s, r> > 01057 * for all rays in the homogeneous space of the two cones that correspond1058 * to the input polyhedra bset1 and bset2.1059 *1060 * We compute s as a vector that satisfies1061 *1062 *	s = \sum_j \alpha_{ij} h_{ij}	for i = 1,2			(*)1063 *1064 * with h_{ij} the normals of the facets of polyhedron i1065 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}1066 * strictly positive numbers.  For simplicity we impose \alpha_{ij} >= 1.1067 * We first set up an LP with as variables the \alpha{ij}.1068 * In this formulation, for each polyhedron i,1069 * the first constraint is the positivity constraint, followed by pairs1070 * of variables for the equalities, followed by variables for the inequalities.1071 * We then simply pick a feasible solution and compute s using (*).1072 *1073 * Note that we simply pick any valid direction and make no attempt1074 * to pick a "good" or even the "best" valid direction.1075 */1076static __isl_give isl_vec *valid_direction(1077	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)1078{1079	struct isl_basic_set *lp;1080	struct isl_tab *tab;1081	struct isl_vec *sample = NULL;1082	struct isl_vec *dir;1083	isl_size d;1084	int i;1085	int n;1086 1087	if (!bset1 || !bset2)1088		goto error;1089	lp = valid_direction_lp(isl_basic_set_copy(bset1),1090				isl_basic_set_copy(bset2));1091	tab = isl_tab_from_basic_set(lp, 0);1092	sample = isl_tab_get_sample_value(tab);1093	isl_tab_free(tab);1094	isl_basic_set_free(lp);1095	if (!sample)1096		goto error;1097	d = isl_basic_set_dim(bset1, isl_dim_all);1098	if (d < 0)1099		goto error;1100	dir = isl_vec_alloc(bset1->ctx, 1 + d);1101	if (!dir)1102		goto error;1103	isl_seq_clr(dir->block.data + 1, dir->size - 1);1104	n = 1;1105	/* positivity constraint 1 >= 0 */1106	isl_int_set(dir->block.data[0], sample->block.data[n]); n++;1107	for (i = 0; i < bset1->n_eq; ++i) {1108		isl_int_sub(sample->block.data[n],1109			    sample->block.data[n], sample->block.data[n+1]);1110		isl_seq_combine(dir->block.data,1111				bset1->ctx->one, dir->block.data,1112				sample->block.data[n], bset1->eq[i], 1 + d);1113 1114		n += 2;1115	}1116	for (i = 0; i < bset1->n_ineq; ++i)1117		isl_seq_combine(dir->block.data,1118				bset1->ctx->one, dir->block.data,1119				sample->block.data[n++], bset1->ineq[i], 1 + d);1120	isl_vec_free(sample);1121	isl_seq_normalize(bset1->ctx, dir->el, dir->size);1122	isl_basic_set_free(bset1);1123	isl_basic_set_free(bset2);1124	return dir;1125error:1126	isl_vec_free(sample);1127	isl_basic_set_free(bset1);1128	isl_basic_set_free(bset2);1129	return NULL;1130}1131 1132/* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},1133 * compute b_i' + A_i' x' >= 0, with1134 *1135 *	[ b_i A_i ]        [ y' ]		              [ y' ]1136 *	[  1   0  ] S^{-1} [ x' ] >= 0	or	[ b_i' A_i' ] [ x' ] >= 01137 *1138 * In particular, add the "positivity constraint" and then perform1139 * the mapping.1140 */1141static __isl_give isl_basic_set *homogeneous_map(__isl_take isl_basic_set *bset,1142	__isl_take isl_mat *T)1143{1144	int k;1145	isl_size total;1146 1147	total = isl_basic_set_dim(bset, isl_dim_all);1148	if (total < 0)1149		goto error;1150	bset = isl_basic_set_extend_constraints(bset, 0, 1);1151	k = isl_basic_set_alloc_inequality(bset);1152	if (k < 0)1153		goto error;1154	isl_seq_clr(bset->ineq[k] + 1, total);1155	isl_int_set_si(bset->ineq[k][0], 1);1156	bset = isl_basic_set_preimage(bset, T);1157	return bset;1158error:1159	isl_mat_free(T);1160	isl_basic_set_free(bset);1161	return NULL;1162}1163 1164/* Compute the convex hull of a pair of basic sets without any parameters or1165 * integer divisions, where the convex hull is known to be pointed,1166 * but the basic sets may be unbounded.1167 *1168 * We turn this problem into the computation of a convex hull of a pair1169 * _bounded_ polyhedra by "changing the direction of the homogeneous1170 * dimension".  This idea is due to Matthias Koeppe.1171 *1172 * Consider the cones in homogeneous space that correspond to the1173 * input polyhedra.  The rays of these cones are also rays of the1174 * polyhedra if the coordinate that corresponds to the homogeneous1175 * dimension is zero.  That is, if the inner product of the rays1176 * with the homogeneous direction is zero.1177 * The cones in the homogeneous space can also be considered to1178 * correspond to other pairs of polyhedra by chosing a different1179 * homogeneous direction.  To ensure that both of these polyhedra1180 * are bounded, we need to make sure that all rays of the cones1181 * correspond to vertices and not to rays.1182 * Let s be a direction such that <s, r> > 0 for all rays r of both cones.1183 * Then using s as a homogeneous direction, we obtain a pair of polytopes.1184 * The vector s is computed in valid_direction.1185 *1186 * Note that we need to consider _all_ rays of the cones and not just1187 * the rays that correspond to rays in the polyhedra.  If we were to1188 * only consider those rays and turn them into vertices, then we1189 * may inadvertently turn some vertices into rays.1190 *1191 * The standard homogeneous direction is the unit vector in the 0th coordinate.1192 * We therefore transform the two polyhedra such that the selected1193 * direction is mapped onto this standard direction and then proceed1194 * with the normal computation.1195 * Let S be a non-singular square matrix with s as its first row,1196 * then we want to map the polyhedra to the space1197 *1198 *	[ y' ]     [ y ]		[ y ]          [ y' ]1199 *	[ x' ] = S [ x ]	i.e.,	[ x ] = S^{-1} [ x' ]1200 *1201 * We take S to be the unimodular completion of s to limit the growth1202 * of the coefficients in the following computations.1203 *1204 * Let b_i + A_i x >= 0 be the constraints of polyhedron i.1205 * We first move to the homogeneous dimension1206 *1207 *	b_i y + A_i x >= 0		[ b_i A_i ] [ y ]    [ 0 ]1208 *	    y         >= 0	or	[  1   0  ] [ x ] >= [ 0 ]1209 *1210 * Then we change directoin1211 *1212 *	[ b_i A_i ]        [ y' ]		              [ y' ]1213 *	[  1   0  ] S^{-1} [ x' ] >= 0	or	[ b_i' A_i' ] [ x' ] >= 01214 *1215 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 01216 * resulting in b' + A' x' >= 0, which we then convert back1217 *1218 *	            [ y ]		        [ y ]1219 *	[ b' A' ] S [ x ] >= 0	or	[ b A ] [ x ] >= 01220 *1221 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.1222 */1223static __isl_give isl_basic_set *convex_hull_pair_pointed(1224	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)1225{1226	struct isl_ctx *ctx = NULL;1227	struct isl_vec *dir = NULL;1228	struct isl_mat *T = NULL;1229	struct isl_mat *T2 = NULL;1230	struct isl_basic_set *hull;1231	struct isl_set *set;1232 1233	if (!bset1 || !bset2)1234		goto error;1235	ctx = isl_basic_set_get_ctx(bset1);1236	dir = valid_direction(isl_basic_set_copy(bset1),1237				isl_basic_set_copy(bset2));1238	if (!dir)1239		goto error;1240	T = isl_mat_alloc(ctx, dir->size, dir->size);1241	if (!T)1242		goto error;1243	isl_seq_cpy(T->row[0], dir->block.data, dir->size);1244	T = isl_mat_unimodular_complete(T, 1);1245	T2 = isl_mat_right_inverse(isl_mat_copy(T));1246 1247	bset1 = homogeneous_map(bset1, isl_mat_copy(T2));1248	bset2 = homogeneous_map(bset2, T2);1249	set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);1250	set = isl_set_add_basic_set(set, bset1);1251	set = isl_set_add_basic_set(set, bset2);1252	hull = uset_convex_hull(set);1253	hull = isl_basic_set_preimage(hull, T);1254	 1255	isl_vec_free(dir);1256 1257	return hull;1258error:1259	isl_vec_free(dir);1260	isl_basic_set_free(bset1);1261	isl_basic_set_free(bset2);1262	return NULL;1263}1264 1265static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set);1266static __isl_give isl_basic_set *modulo_affine_hull(1267	__isl_take isl_set *set, __isl_take isl_basic_set *affine_hull);1268 1269/* Compute the convex hull of a pair of basic sets without any parameters or1270 * integer divisions.1271 *1272 * This function is called from uset_convex_hull_unbounded, which1273 * means that the complete convex hull is unbounded.  Some pairs1274 * of basic sets may still be bounded, though.1275 * They may even lie inside a lower dimensional space, in which1276 * case they need to be handled inside their affine hull since1277 * the main algorithm assumes that the result is full-dimensional.1278 *1279 * If the convex hull of the two basic sets would have a non-trivial1280 * lineality space, we first project out this lineality space.1281 */1282static __isl_give isl_basic_set *convex_hull_pair(1283	__isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)1284{1285	isl_basic_set *lin, *aff;1286	isl_bool bounded1, bounded2;1287	isl_size total;1288 1289	if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM)1290		return convex_hull_pair_elim(bset1, bset2);1291 1292	aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),1293						    isl_basic_set_copy(bset2)));1294	if (!aff)1295		goto error;1296	if (aff->n_eq != 0) 1297		return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);1298	isl_basic_set_free(aff);1299 1300	bounded1 = isl_basic_set_is_bounded(bset1);1301	bounded2 = isl_basic_set_is_bounded(bset2);1302 1303	if (bounded1 < 0 || bounded2 < 0)1304		goto error;1305 1306	if (bounded1 && bounded2)1307		return uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));1308 1309	if (bounded1 || bounded2)1310		return convex_hull_pair_pointed(bset1, bset2);1311 1312	lin = induced_lineality_space(isl_basic_set_copy(bset1),1313				      isl_basic_set_copy(bset2));1314	if (!lin)1315		goto error;1316	if (isl_basic_set_plain_is_universe(lin)) {1317		isl_basic_set_free(bset1);1318		isl_basic_set_free(bset2);1319		return lin;1320	}1321	total = isl_basic_set_dim(lin, isl_dim_all);1322	if (lin->n_eq < total) {1323		struct isl_set *set;1324		set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);1325		set = isl_set_add_basic_set(set, bset1);1326		set = isl_set_add_basic_set(set, bset2);1327		return modulo_lineality(set, lin);1328	}1329	isl_basic_set_free(lin);1330	if (total < 0)1331		goto error;1332 1333	return convex_hull_pair_pointed(bset1, bset2);1334error:1335	isl_basic_set_free(bset1);1336	isl_basic_set_free(bset2);1337	return NULL;1338}1339 1340/* Compute the lineality space of a basic set.1341 * We basically just drop the constants and turn every inequality1342 * into an equality.1343 * Any explicit representations of local variables are removed1344 * because they may no longer be valid representations1345 * in the lineality space.1346 */1347__isl_give isl_basic_set *isl_basic_set_lineality_space(1348	__isl_take isl_basic_set *bset)1349{1350	int i, k;1351	struct isl_basic_set *lin = NULL;1352	isl_size n_div, dim;1353 1354	n_div = isl_basic_set_dim(bset, isl_dim_div);1355	dim = isl_basic_set_dim(bset, isl_dim_all);1356	if (n_div < 0 || dim < 0)1357		return isl_basic_set_free(bset);1358 1359	lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset),1360					n_div, dim, 0);1361	for (i = 0; i < n_div; ++i)1362		if (isl_basic_set_alloc_div(lin) < 0)1363			goto error;1364	if (!lin)1365		goto error;1366	for (i = 0; i < bset->n_eq; ++i) {1367		k = isl_basic_set_alloc_equality(lin);1368		if (k < 0)1369			goto error;1370		isl_int_set_si(lin->eq[k][0], 0);1371		isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim);1372	}1373	lin = isl_basic_set_gauss(lin, NULL);1374	if (!lin)1375		goto error;1376	for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) {1377		k = isl_basic_set_alloc_equality(lin);1378		if (k < 0)1379			goto error;1380		isl_int_set_si(lin->eq[k][0], 0);1381		isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim);1382		lin = isl_basic_set_gauss(lin, NULL);1383		if (!lin)1384			goto error;1385	}1386	isl_basic_set_free(bset);1387	return lin;1388error:1389	isl_basic_set_free(lin);1390	isl_basic_set_free(bset);1391	return NULL;1392}1393 1394/* Compute the (linear) hull of the lineality spaces of the basic sets in the1395 * set "set".1396 */1397__isl_give isl_basic_set *isl_set_combined_lineality_space(1398	__isl_take isl_set *set)1399{1400	int i;1401	struct isl_set *lin = NULL;1402 1403	if (!set)1404		return NULL;1405	if (set->n == 0) {1406		isl_space *space = isl_set_get_space(set);1407		isl_set_free(set);1408		return isl_basic_set_empty(space);1409	}1410 1411	lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0);1412	for (i = 0; i < set->n; ++i)1413		lin = isl_set_add_basic_set(lin,1414		    isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));1415	isl_set_free(set);1416	return isl_set_affine_hull(lin);1417}1418 1419/* Compute the convex hull of a set without any parameters or1420 * integer divisions.1421 * In each step, we combined two basic sets until only one1422 * basic set is left.1423 * The input basic sets are assumed not to have a non-trivial1424 * lineality space.  If any of the intermediate results has1425 * a non-trivial lineality space, it is projected out.1426 */1427static __isl_give isl_basic_set *uset_convex_hull_unbounded(1428	__isl_take isl_set *set)1429{1430	isl_basic_set_list *list;1431 1432	list = isl_set_get_basic_set_list(set);1433	isl_set_free(set);1434 1435	while (list) {1436		isl_size n, total;1437		struct isl_basic_set *t;1438		isl_basic_set *bset1, *bset2;1439 1440		n = isl_basic_set_list_n_basic_set(list);1441		if (n < 0)1442			goto error;1443		if (n < 2)1444			isl_die(isl_basic_set_list_get_ctx(list),1445				isl_error_internal,1446				"expecting at least two elements", goto error);1447		bset1 = isl_basic_set_list_get_basic_set(list, n - 1);1448		bset2 = isl_basic_set_list_get_basic_set(list, n - 2);1449		bset1 = convex_hull_pair(bset1, bset2);1450		if (n == 2) {1451			isl_basic_set_list_free(list);1452			return bset1;1453		}1454		bset1 = isl_basic_set_underlying_set(bset1);1455		list = isl_basic_set_list_drop(list, n - 2, 2);1456		list = isl_basic_set_list_add(list, bset1);1457 1458		t = isl_basic_set_list_get_basic_set(list, n - 2);1459		t = isl_basic_set_lineality_space(t);1460		if (!t)1461			goto error;1462		if (isl_basic_set_plain_is_universe(t)) {1463			isl_basic_set_list_free(list);1464			return t;1465		}1466		total = isl_basic_set_dim(t, isl_dim_all);1467		if (t->n_eq < total) {1468			set = isl_basic_set_list_union(list);1469			return modulo_lineality(set, t);1470		}1471		isl_basic_set_free(t);1472		if (total < 0)1473			goto error;1474	}1475 1476	return NULL;1477error:1478	isl_basic_set_list_free(list);1479	return NULL;1480}1481 1482/* Compute an initial hull for wrapping containing a single initial1483 * facet.1484 * This function assumes that the given set is bounded.1485 */1486static __isl_give isl_basic_set *initial_hull(__isl_take isl_basic_set *hull,1487	__isl_keep isl_set *set)1488{1489	struct isl_mat *bounds = NULL;1490	isl_size dim;1491	int k;1492 1493	if (!hull)1494		goto error;1495	bounds = initial_facet_constraint(set);1496	if (!bounds)1497		goto error;1498	k = isl_basic_set_alloc_inequality(hull);1499	if (k < 0)1500		goto error;1501	dim = isl_set_dim(set, isl_dim_set);1502	if (dim < 0)1503		goto error;1504	isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);1505	isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col);1506	isl_mat_free(bounds);1507 1508	return hull;1509error:1510	isl_basic_set_free(hull);1511	isl_mat_free(bounds);1512	return NULL;1513}1514 1515struct max_constraint {1516	struct isl_mat *c;1517	int	 	count;1518	int		ineq;1519};1520 1521static isl_bool max_constraint_equal(const void *entry, const void *val)1522{1523	struct max_constraint *a = (struct max_constraint *)entry;1524	isl_int *b = (isl_int *)val;1525 1526	return isl_bool_ok(isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1));1527}1528 1529static isl_stat update_constraint(struct isl_ctx *ctx,1530	struct isl_hash_table *table,1531	isl_int *con, unsigned len, int n, int ineq)1532{1533	struct isl_hash_table_entry *entry;1534	struct max_constraint *c;1535	uint32_t c_hash;1536 1537	c_hash = isl_seq_get_hash(con + 1, len);1538	entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,1539			con + 1, 0);1540	if (!entry)1541		return isl_stat_error;1542	if (entry == isl_hash_table_entry_none)1543		return isl_stat_ok;1544	c = entry->data;1545	if (c->count < n) {1546		isl_hash_table_remove(ctx, table, entry);1547		return isl_stat_ok;1548	}1549	c->count++;1550	if (isl_int_gt(c->c->row[0][0], con[0]))1551		return isl_stat_ok;1552	if (isl_int_eq(c->c->row[0][0], con[0])) {1553		if (ineq)1554			c->ineq = ineq;1555		return isl_stat_ok;1556	}1557	c->c = isl_mat_cow(c->c);1558	isl_int_set(c->c->row[0][0], con[0]);1559	c->ineq = ineq;1560 1561	return isl_stat_ok;1562}1563 1564/* Check whether the constraint hash table "table" contains the constraint1565 * "con".1566 */1567static isl_bool has_constraint(struct isl_ctx *ctx,1568	struct isl_hash_table *table, isl_int *con, unsigned len, int n)1569{1570	struct isl_hash_table_entry *entry;1571	struct max_constraint *c;1572	uint32_t c_hash;1573 1574	c_hash = isl_seq_get_hash(con + 1, len);1575	entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,1576			con + 1, 0);1577	if (!entry)1578		return isl_bool_error;1579	if (entry == isl_hash_table_entry_none)1580		return isl_bool_false;1581	c = entry->data;1582	if (c->count < n)1583		return isl_bool_false;1584	return isl_bool_ok(isl_int_eq(c->c->row[0][0], con[0]));1585}1586 1587/* Are the constraints of "bset" known to be facets?1588 * If there are any equality constraints, then they are not.1589 * If there may be redundant constraints, then those1590 * redundant constraints are not facets.1591 */1592static isl_bool has_facets(__isl_keep isl_basic_set *bset)1593{1594	isl_size n_eq;1595 1596	n_eq = isl_basic_set_n_equality(bset);1597	if (n_eq < 0)1598		return isl_bool_error;1599	if (n_eq != 0)1600		return isl_bool_false;1601	return ISL_F_ISSET(bset, ISL_BASIC_SET_NO_REDUNDANT);1602}1603 1604/* Check for inequality constraints of a basic set without equalities1605 * or redundant constraints1606 * such that the same or more stringent copies of the constraint appear1607 * in all of the basic sets.  Such constraints are necessarily facet1608 * constraints of the convex hull.1609 *1610 * If the resulting basic set is by chance identical to one of1611 * the basic sets in "set", then we know that this basic set contains1612 * all other basic sets and is therefore the convex hull of set.1613 * In this case we set *is_hull to 1.1614 */1615static __isl_give isl_basic_set *common_constraints(1616	__isl_take isl_basic_set *hull, __isl_keep isl_set *set, int *is_hull)1617{1618	int i, j, s, n;1619	int min_constraints;1620	int best;1621	struct max_constraint *constraints = NULL;1622	struct isl_hash_table *table = NULL;1623	isl_size total;1624 1625	*is_hull = 0;1626 1627	for (i = 0; i < set->n; ++i) {1628		isl_bool facets = has_facets(set->p[i]);1629		if (facets < 0)1630			return isl_basic_set_free(hull);1631		if (facets)1632			break;1633	}1634	if (i >= set->n)1635		return hull;1636	min_constraints = set->p[i]->n_ineq;1637	best = i;1638	for (i = best + 1; i < set->n; ++i) {1639		isl_bool facets = has_facets(set->p[i]);1640		if (facets < 0)1641			return isl_basic_set_free(hull);1642		if (!facets)1643			continue;1644		if (set->p[i]->n_ineq >= min_constraints)1645			continue;1646		min_constraints = set->p[i]->n_ineq;1647		best = i;1648	}1649	constraints = isl_calloc_array(hull->ctx, struct max_constraint,1650					min_constraints);1651	if (!constraints)1652		return hull;1653	table = isl_alloc_type(hull->ctx, struct isl_hash_table);1654	if (isl_hash_table_init(hull->ctx, table, min_constraints))1655		goto error;1656 1657	total = isl_set_dim(set, isl_dim_all);1658	if (total < 0)1659		goto error;1660	for (i = 0; i < set->p[best]->n_ineq; ++i) {1661		constraints[i].c = isl_mat_sub_alloc6(hull->ctx,1662			set->p[best]->ineq + i, 0, 1, 0, 1 + total);1663		if (!constraints[i].c)1664			goto error;1665		constraints[i].ineq = 1;1666	}1667	for (i = 0; i < min_constraints; ++i) {1668		struct isl_hash_table_entry *entry;1669		uint32_t c_hash;1670		c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);1671		entry = isl_hash_table_find(hull->ctx, table, c_hash,1672			max_constraint_equal, constraints[i].c->row[0] + 1, 1);1673		if (!entry)1674			goto error;1675		isl_assert(hull->ctx, !entry->data, goto error);1676		entry->data = &constraints[i];1677	}1678 1679	n = 0;1680	for (s = 0; s < set->n; ++s) {1681		if (s == best)1682			continue;1683 1684		for (i = 0; i < set->p[s]->n_eq; ++i) {1685			isl_int *eq = set->p[s]->eq[i];1686			for (j = 0; j < 2; ++j) {1687				isl_seq_neg(eq, eq, 1 + total);1688				if (update_constraint(hull->ctx, table,1689						    eq, total, n, 0) < 0)1690					goto error;1691			}1692		}1693		for (i = 0; i < set->p[s]->n_ineq; ++i) {1694			isl_int *ineq = set->p[s]->ineq[i];1695			if (update_constraint(hull->ctx, table, ineq, total, n,1696					set->p[s]->n_eq == 0) < 0)1697				goto error;1698		}1699		++n;1700	}1701 1702	for (i = 0; i < min_constraints; ++i) {1703		if (constraints[i].count < n)1704			continue;1705		if (!constraints[i].ineq)1706			continue;1707		j = isl_basic_set_alloc_inequality(hull);1708		if (j < 0)1709			goto error;1710		isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total);1711	}1712 1713	for (s = 0; s < set->n; ++s) {1714		if (set->p[s]->n_eq)1715			continue;1716		if (set->p[s]->n_ineq != hull->n_ineq)1717			continue;1718		for (i = 0; i < set->p[s]->n_ineq; ++i) {1719			isl_bool has;1720			isl_int *ineq = set->p[s]->ineq[i];1721			has = has_constraint(hull->ctx, table, ineq, total, n);1722			if (has < 0)1723				goto error;1724			if (!has)1725				break;1726		}1727		if (i == set->p[s]->n_ineq)1728			*is_hull = 1;1729	}1730 1731	isl_hash_table_clear(table);1732	for (i = 0; i < min_constraints; ++i)1733		isl_mat_free(constraints[i].c);1734	free(constraints);1735	free(table);1736	return hull;1737error:1738	isl_hash_table_clear(table);1739	free(table);1740	if (constraints)1741		for (i = 0; i < min_constraints; ++i)1742			isl_mat_free(constraints[i].c);1743	free(constraints);1744	return hull;1745}1746 1747/* Create a template for the convex hull of "set" and fill it up1748 * obvious facet constraints, if any.  If the result happens to1749 * be the convex hull of "set" then *is_hull is set to 1.1750 */1751static __isl_give isl_basic_set *proto_hull(__isl_keep isl_set *set,1752	int *is_hull)1753{1754	struct isl_basic_set *hull;1755	unsigned n_ineq;1756	int i;1757 1758	n_ineq = 1;1759	for (i = 0; i < set->n; ++i) {1760		n_ineq += set->p[i]->n_eq;1761		n_ineq += set->p[i]->n_ineq;1762	}1763	hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);1764	hull = isl_basic_set_set_rational(hull);1765	if (!hull)1766		return NULL;1767	return common_constraints(hull, set, is_hull);1768}1769 1770static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set)1771{1772	struct isl_basic_set *hull;1773	int is_hull;1774 1775	hull = proto_hull(set, &is_hull);1776	if (hull && !is_hull) {1777		if (hull->n_ineq == 0)1778			hull = initial_hull(hull, set);1779		hull = extend(hull, set);1780	}1781	isl_set_free(set);1782 1783	return hull;1784}1785 1786/* Compute the convex hull of a set without any parameters or1787 * integer divisions.  Depending on whether the set is bounded,1788 * we pass control to the wrapping based convex hull or1789 * the Fourier-Motzkin elimination based convex hull.1790 * We also handle a few special cases before checking the boundedness.1791 */1792static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set)1793{1794	isl_bool bounded;1795	isl_size dim;1796	struct isl_basic_set *convex_hull = NULL;1797	struct isl_basic_set *lin;1798 1799	dim = isl_set_dim(set, isl_dim_all);1800	if (dim < 0)1801		goto error;1802	if (dim == 0)1803		return convex_hull_0d(set);1804 1805	set = isl_set_coalesce(set);1806	set = isl_set_set_rational(set);1807 1808	if (!set)1809		return NULL;1810	if (set->n == 1) {1811		convex_hull = isl_basic_set_copy(set->p[0]);1812		isl_set_free(set);1813		return convex_hull;1814	}1815	if (dim == 1)1816		return convex_hull_1d(set);1817 1818	bounded = isl_set_is_bounded(set);1819	if (bounded < 0)1820		goto error;1821	if (bounded && set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP)1822		return uset_convex_hull_wrap(set);1823 1824	lin = isl_set_combined_lineality_space(isl_set_copy(set));1825	if (!lin)1826		goto error;1827	if (isl_basic_set_plain_is_universe(lin)) {1828		isl_set_free(set);1829		return lin;1830	}1831	if (lin->n_eq < dim)1832		return modulo_lineality(set, lin);1833	isl_basic_set_free(lin);1834 1835	return uset_convex_hull_unbounded(set);1836error:1837	isl_set_free(set);1838	isl_basic_set_free(convex_hull);1839	return NULL;1840}1841 1842/* This is the core procedure, where "set" is a "pure" set, i.e.,1843 * without parameters or divs and where the convex hull of set is1844 * known to be full-dimensional.1845 */1846static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(1847	__isl_take isl_set *set)1848{1849	struct isl_basic_set *convex_hull = NULL;1850	isl_size dim;1851 1852	dim = isl_set_dim(set, isl_dim_all);1853	if (dim < 0)1854		goto error;1855 1856	if (dim == 0) {1857		convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));1858		isl_set_free(set);1859		convex_hull = isl_basic_set_set_rational(convex_hull);1860		return convex_hull;1861	}1862 1863	set = isl_set_set_rational(set);1864	set = isl_set_coalesce(set);1865	if (!set)1866		goto error;1867	if (set->n == 1) {1868		convex_hull = isl_basic_set_copy(set->p[0]);1869		isl_set_free(set);1870		convex_hull = isl_basic_map_remove_redundancies(convex_hull);1871		return convex_hull;1872	}1873	if (dim == 1)1874		return convex_hull_1d(set);1875 1876	return uset_convex_hull_wrap(set);1877error:1878	isl_set_free(set);1879	return NULL;1880}1881 1882/* Compute the convex hull of set "set" with affine hull "affine_hull",1883 * We first remove the equalities (transforming the set), compute the1884 * convex hull of the transformed set and then add the equalities back1885 * (after performing the inverse transformation.1886 */1887static __isl_give isl_basic_set *modulo_affine_hull(1888	__isl_take isl_set *set, __isl_take isl_basic_set *affine_hull)1889{1890	struct isl_mat *T;1891	struct isl_mat *T2;1892	struct isl_basic_set *dummy;1893	struct isl_basic_set *convex_hull;1894 1895	dummy = isl_basic_set_remove_equalities(1896			isl_basic_set_copy(affine_hull), &T, &T2);1897	if (!dummy)1898		goto error;1899	isl_basic_set_free(dummy);1900	set = isl_set_preimage(set, T);1901	convex_hull = uset_convex_hull(set);1902	convex_hull = isl_basic_set_preimage(convex_hull, T2);1903	convex_hull = isl_basic_set_intersect(convex_hull, affine_hull);1904	return convex_hull;1905error:1906	isl_mat_free(T);1907	isl_mat_free(T2);1908	isl_basic_set_free(affine_hull);1909	isl_set_free(set);1910	return NULL;1911}1912 1913/* Return an empty basic map living in the same space as "map".1914 */1915static __isl_give isl_basic_map *replace_map_by_empty_basic_map(1916	__isl_take isl_map *map)1917{1918	isl_space *space;1919 1920	space = isl_map_get_space(map);1921	isl_map_free(map);1922	return isl_basic_map_empty(space);1923}1924 1925/* Compute the convex hull of a map.1926 *1927 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,1928 * specifically, the wrapping of facets to obtain new facets.1929 */1930__isl_give isl_basic_map *isl_map_convex_hull(__isl_take isl_map *map)1931{1932	struct isl_basic_set *bset;1933	struct isl_basic_map *model = NULL;1934	struct isl_basic_set *affine_hull = NULL;1935	struct isl_basic_map *convex_hull = NULL;1936	struct isl_set *set = NULL;1937 1938	map = isl_map_detect_equalities(map);1939	map = isl_map_align_divs_internal(map);1940	if (!map)1941		goto error;1942 1943	if (map->n == 0)1944		return replace_map_by_empty_basic_map(map);1945 1946	model = isl_basic_map_copy(map->p[0]);1947	set = isl_map_underlying_set(map);1948	if (!set)1949		goto error;1950 1951	affine_hull = isl_set_affine_hull(isl_set_copy(set));1952	if (!affine_hull)1953		goto error;1954	if (affine_hull->n_eq != 0)1955		bset = modulo_affine_hull(set, affine_hull);1956	else {1957		isl_basic_set_free(affine_hull);1958		bset = uset_convex_hull(set);1959	}1960 1961	convex_hull = isl_basic_map_overlying_set(bset, model);1962	if (!convex_hull)1963		return NULL;1964 1965	ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);1966	ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);1967	ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);1968	return convex_hull;1969error:1970	isl_set_free(set);1971	isl_basic_map_free(model);1972	return NULL;1973}1974 1975__isl_give isl_basic_set *isl_set_convex_hull(__isl_take isl_set *set)1976{1977	return bset_from_bmap(isl_map_convex_hull(set_to_map(set)));1978}1979 1980__isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)1981{1982	isl_basic_map *hull;1983 1984	hull = isl_map_convex_hull(map);1985	return isl_basic_map_remove_divs(hull);1986}1987 1988__isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)1989{1990	return bset_from_bmap(isl_map_polyhedral_hull(set_to_map(set)));1991}1992 1993struct sh_data_entry {1994	struct isl_hash_table	*table;1995	struct isl_tab		*tab;1996};1997 1998/* Holds the data needed during the simple hull computation.1999 * In particular,2000 *	n		the number of basic sets in the original set2001 *	hull_table	a hash table of already computed constraints2002 *			in the simple hull2003 *	p		for each basic set,2004 *		table		a hash table of the constraints2005 *		tab		the tableau corresponding to the basic set2006 */2007struct sh_data {2008	struct isl_ctx		*ctx;2009	unsigned		n;2010	struct isl_hash_table	*hull_table;2011	struct sh_data_entry	p[1];2012};2013 2014static void sh_data_free(struct sh_data *data)2015{2016	int i;2017 2018	if (!data)2019		return;2020	isl_hash_table_free(data->ctx, data->hull_table);2021	for (i = 0; i < data->n; ++i) {2022		isl_hash_table_free(data->ctx, data->p[i].table);2023		isl_tab_free(data->p[i].tab);2024	}2025	free(data);2026}2027 2028struct ineq_cmp_data {2029	unsigned	len;2030	isl_int		*p;2031};2032 2033static isl_bool has_ineq(const void *entry, const void *val)2034{2035	isl_int *row = (isl_int *)entry;2036	struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;2037 2038	return isl_bool_ok(isl_seq_eq(row + 1, v->p + 1, v->len) ||2039			   isl_seq_is_neg(row + 1, v->p + 1, v->len));2040}2041 2042static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,2043			isl_int *ineq, unsigned len)2044{2045	uint32_t c_hash;2046	struct ineq_cmp_data v;2047	struct isl_hash_table_entry *entry;2048 2049	v.len = len;2050	v.p = ineq;2051	c_hash = isl_seq_get_hash(ineq + 1, len);2052	entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);2053	if (!entry)2054		return - 1;2055	entry->data = ineq;2056	return 0;2057}2058 2059/* Fill hash table "table" with the constraints of "bset".2060 * Equalities are added as two inequalities.2061 * The value in the hash table is a pointer to the (in)equality of "bset".2062 */2063static isl_stat hash_basic_set(struct isl_hash_table *table,2064	__isl_keep isl_basic_set *bset)2065{2066	int i, j;2067	isl_size dim = isl_basic_set_dim(bset, isl_dim_all);2068 2069	if (dim < 0)2070		return isl_stat_error;2071	for (i = 0; i < bset->n_eq; ++i) {2072		for (j = 0; j < 2; ++j) {2073			isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim);2074			if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0)2075				return isl_stat_error;2076		}2077	}2078	for (i = 0; i < bset->n_ineq; ++i) {2079		if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0)2080			return isl_stat_error;2081	}2082	return isl_stat_ok;2083}2084 2085static struct sh_data *sh_data_alloc(__isl_keep isl_set *set, unsigned n_ineq)2086{2087	struct sh_data *data;2088	int i;2089 2090	data = isl_calloc(set->ctx, struct sh_data,2091		sizeof(struct sh_data) +2092		(set->n - 1) * sizeof(struct sh_data_entry));2093	if (!data)2094		return NULL;2095	data->ctx = set->ctx;2096	data->n = set->n;2097	data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq);2098	if (!data->hull_table)2099		goto error;2100	for (i = 0; i < set->n; ++i) {2101		data->p[i].table = isl_hash_table_alloc(set->ctx,2102				    2 * set->p[i]->n_eq + set->p[i]->n_ineq);2103		if (!data->p[i].table)2104			goto error;2105		if (hash_basic_set(data->p[i].table, set->p[i]) < 0)2106			goto error;2107	}2108	return data;2109error:2110	sh_data_free(data);2111	return NULL;2112}2113 2114/* Check if inequality "ineq" is a bound for basic set "j" or if2115 * it can be relaxed (by increasing the constant term) to become2116 * a bound for that basic set.  In the latter case, the constant2117 * term is updated.2118 * Relaxation of the constant term is only allowed if "shift" is set.2119 *2120 * Return 1 if "ineq" is a bound2121 *	  0 if "ineq" may attain arbitrarily small values on basic set "j"2122 *	 -1 if some error occurred2123 */2124static int is_bound(struct sh_data *data, __isl_keep isl_set *set, int j,2125	isl_int *ineq, int shift)2126{2127	enum isl_lp_result res;2128	isl_int opt;2129 2130	if (!data->p[j].tab) {2131		data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0);2132		if (!data->p[j].tab)2133			return -1;2134	}2135 2136	isl_int_init(opt);2137 2138	res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one,2139				&opt, NULL, 0);2140	if (res == isl_lp_ok && isl_int_is_neg(opt)) {2141		if (shift)2142			isl_int_sub(ineq[0], ineq[0], opt);2143		else2144			res = isl_lp_unbounded;2145	}2146 2147	isl_int_clear(opt);2148 2149	return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :2150	       res == isl_lp_unbounded ? 0 : -1;2151}2152 2153/* Set the constant term of "ineq" to the maximum of those of the constraints2154 * in the basic sets of "set" following "i" that are parallel to "ineq".2155 * That is, if any of the basic sets of "set" following "i" have a more2156 * relaxed copy of "ineq", then replace "ineq" by the most relaxed copy.2157 * "c_hash" is the hash value of the linear part of "ineq".2158 * "v" has been set up for use by has_ineq.2159 *2160 * Note that the two inequality constraints corresponding to an equality are2161 * represented by the same inequality constraint in data->p[j].table2162 * (but with different hash values).  This means the constraint (or at2163 * least its constant term) may need to be temporarily negated to get2164 * the actually hashed constraint.2165 */2166static isl_stat set_max_constant_term(struct sh_data *data,2167	__isl_keep isl_set *set,2168	int i, isl_int *ineq, uint32_t c_hash, struct ineq_cmp_data *v)2169{2170	int j;2171	isl_ctx *ctx;2172	struct isl_hash_table_entry *entry;2173 2174	ctx = isl_set_get_ctx(set);2175	for (j = i + 1; j < set->n; ++j) {2176		int neg;2177		isl_int *ineq_j;2178 2179		entry = isl_hash_table_find(ctx, data->p[j].table,2180						c_hash, &has_ineq, v, 0);2181		if (!entry)2182			return isl_stat_error;2183		if (entry == isl_hash_table_entry_none)2184			continue;2185 2186		ineq_j = entry->data;2187		neg = isl_seq_is_neg(ineq_j + 1, ineq + 1, v->len);2188		if (neg)2189			isl_int_neg(ineq_j[0], ineq_j[0]);2190		if (isl_int_gt(ineq_j[0], ineq[0]))2191			isl_int_set(ineq[0], ineq_j[0]);2192		if (neg)2193			isl_int_neg(ineq_j[0], ineq_j[0]);2194	}2195 2196	return isl_stat_ok;2197}2198 2199/* Check if inequality "ineq" from basic set "i" is or can be relaxed to2200 * become a bound on the whole set.  If so, add the (relaxed) inequality2201 * to "hull".  Relaxation is only allowed if "shift" is set.2202 *2203 * We first check if "hull" already contains a translate of the inequality.2204 * If so, we are done.2205 * Then, we check if any of the previous basic sets contains a translate2206 * of the inequality.  If so, then we have already considered this2207 * inequality and we are done.2208 * Otherwise, for each basic set other than "i", we check if the inequality2209 * is a bound on the basic set, but first replace the constant term2210 * by the maximal value of any translate of the inequality in any2211 * of the following basic sets.2212 * For previous basic sets, we know that they do not contain a translate2213 * of the inequality, so we directly call is_bound.2214 * For following basic sets, we first check if a translate of the2215 * inequality appears in its description.  If so, the constant term2216 * of the inequality has already been updated with respect to this2217 * translate and the inequality is therefore known to be a bound2218 * of this basic set.2219 */2220static __isl_give isl_basic_set *add_bound(__isl_take isl_basic_set *hull,2221	struct sh_data *data, __isl_keep isl_set *set, int i, isl_int *ineq,2222	int shift)2223{2224	uint32_t c_hash;2225	struct ineq_cmp_data v;2226	struct isl_hash_table_entry *entry;2227	int j, k;2228	isl_size total;2229 2230	total = isl_basic_set_dim(hull, isl_dim_all);2231	if (total < 0)2232		return isl_basic_set_free(hull);2233 2234	v.len = total;2235	v.p = ineq;2236	c_hash = isl_seq_get_hash(ineq + 1, v.len);2237 2238	entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,2239					has_ineq, &v, 0);2240	if (!entry)2241		return isl_basic_set_free(hull);2242	if (entry != isl_hash_table_entry_none)2243		return hull;2244 2245	for (j = 0; j < i; ++j) {2246		entry = isl_hash_table_find(hull->ctx, data->p[j].table,2247						c_hash, has_ineq, &v, 0);2248		if (!entry)2249			return isl_basic_set_free(hull);2250		if (entry != isl_hash_table_entry_none)2251			break;2252	}2253	if (j < i)2254		return hull;2255 2256	k = isl_basic_set_alloc_inequality(hull);2257	if (k < 0)2258		goto error;2259	isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);2260 2261	if (set_max_constant_term(data, set, i, hull->ineq[k], c_hash, &v) < 0)2262		goto error;2263	for (j = 0; j < i; ++j) {2264		int bound;2265		bound = is_bound(data, set, j, hull->ineq[k], shift);2266		if (bound < 0)2267			goto error;2268		if (!bound)2269			break;2270	}2271	if (j < i)2272		return isl_basic_set_free_inequality(hull, 1);2273 2274	for (j = i + 1; j < set->n; ++j) {2275		int bound;2276		entry = isl_hash_table_find(hull->ctx, data->p[j].table,2277						c_hash, has_ineq, &v, 0);2278		if (!entry)2279			return isl_basic_set_free(hull);2280		if (entry != isl_hash_table_entry_none)2281			continue;2282		bound = is_bound(data, set, j, hull->ineq[k], shift);2283		if (bound < 0)2284			goto error;2285		if (!bound)2286			break;2287	}2288	if (j < set->n)2289		return isl_basic_set_free_inequality(hull, 1);2290 2291	entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,2292					has_ineq, &v, 1);2293	if (!entry)2294		goto error;2295	entry->data = hull->ineq[k];2296 2297	return hull;2298error:2299	isl_basic_set_free(hull);2300	return NULL;2301}2302 2303/* Check if any inequality from basic set "i" is or can be relaxed to2304 * become a bound on the whole set.  If so, add the (relaxed) inequality2305 * to "hull".  Relaxation is only allowed if "shift" is set.2306 */2307static __isl_give isl_basic_set *add_bounds(__isl_take isl_basic_set *bset,2308	struct sh_data *data, __isl_keep isl_set *set, int i, int shift)2309{2310	int j, k;2311	isl_size dim = isl_basic_set_dim(bset, isl_dim_all);2312 2313	if (dim < 0)2314		return isl_basic_set_free(bset);2315 2316	for (j = 0; j < set->p[i]->n_eq; ++j) {2317		for (k = 0; k < 2; ++k) {2318			isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);2319			bset = add_bound(bset, data, set, i, set->p[i]->eq[j],2320					    shift);2321		}2322	}2323	for (j = 0; j < set->p[i]->n_ineq; ++j)2324		bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift);2325	return bset;2326}2327 2328/* Compute a superset of the convex hull of set that is described2329 * by only (translates of) the constraints in the constituents of set.2330 * Translation is only allowed if "shift" is set.2331 */2332static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,2333	int shift)2334{2335	struct sh_data *data = NULL;2336	struct isl_basic_set *hull = NULL;2337	unsigned n_ineq;2338	int i;2339 2340	if (!set)2341		return NULL;2342 2343	n_ineq = 0;2344	for (i = 0; i < set->n; ++i) {2345		if (!set->p[i])2346			goto error;2347		n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;2348	}2349 2350	hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);2351	if (!hull)2352		goto error;2353 2354	data = sh_data_alloc(set, n_ineq);2355	if (!data)2356		goto error;2357 2358	for (i = 0; i < set->n; ++i)2359		hull = add_bounds(hull, data, set, i, shift);2360 2361	sh_data_free(data);2362	isl_set_free(set);2363 2364	return hull;2365error:2366	sh_data_free(data);2367	isl_basic_set_free(hull);2368	isl_set_free(set);2369	return NULL;2370}2371 2372/* Compute a superset of the convex hull of map that is described2373 * by only (translates of) the constraints in the constituents of map.2374 * Handle trivial cases where map is NULL or contains at most one disjunct.2375 */2376static __isl_give isl_basic_map *map_simple_hull_trivial(2377	__isl_take isl_map *map)2378{2379	isl_basic_map *hull;2380 2381	if (!map)2382		return NULL;2383	if (map->n == 0)2384		return replace_map_by_empty_basic_map(map);2385 2386	hull = isl_basic_map_copy(map->p[0]);2387	isl_map_free(map);2388	return hull;2389}2390 2391/* Return a copy of the simple hull cached inside "map".2392 * "shift" determines whether to return the cached unshifted or shifted2393 * simple hull.2394 */2395static __isl_give isl_basic_map *cached_simple_hull(__isl_take isl_map *map,2396	int shift)2397{2398	isl_basic_map *hull;2399 2400	hull = isl_basic_map_copy(map->cached_simple_hull[shift]);2401	isl_map_free(map);2402 2403	return hull;2404}2405 2406/* Compute a superset of the convex hull of map that is described2407 * by only (translates of) the constraints in the constituents of map.2408 * Translation is only allowed if "shift" is set.2409 *2410 * The constraints are sorted while removing redundant constraints2411 * in order to indicate a preference of which constraints should2412 * be preserved.  In particular, pairs of constraints that are2413 * sorted together are preferred to either both be preserved2414 * or both be removed.  The sorting is performed inside2415 * isl_basic_map_remove_redundancies.2416 *2417 * The result of the computation is stored in map->cached_simple_hull[shift]2418 * such that it can be reused in subsequent calls.  The cache is cleared2419 * whenever the map is modified (in isl_map_cow).2420 * Note that the results need to be stored in the input map for there2421 * to be any chance that they may get reused.  In particular, they2422 * are stored in a copy of the input map that is saved before2423 * the integer division alignment.2424 */2425static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,2426	int shift)2427{2428	struct isl_set *set = NULL;2429	struct isl_basic_map *model = NULL;2430	struct isl_basic_map *hull;2431	struct isl_basic_map *affine_hull;2432	struct isl_basic_set *bset = NULL;2433	isl_map *input;2434 2435	if (!map || map->n <= 1)2436		return map_simple_hull_trivial(map);2437 2438	if (map->cached_simple_hull[shift])2439		return cached_simple_hull(map, shift);2440 2441	map = isl_map_detect_equalities(map);2442	if (!map || map->n <= 1)2443		return map_simple_hull_trivial(map);2444	affine_hull = isl_map_affine_hull(isl_map_copy(map));2445	input = isl_map_copy(map);2446	map = isl_map_align_divs_internal(map);2447	model = map ? isl_basic_map_copy(map->p[0]) : NULL;2448 2449	set = isl_map_underlying_set(map);2450 2451	bset = uset_simple_hull(set, shift);2452 2453	hull = isl_basic_map_overlying_set(bset, model);2454 2455	hull = isl_basic_map_intersect(hull, affine_hull);2456	hull = isl_basic_map_remove_redundancies(hull);2457 2458	if (hull) {2459		ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);2460		ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);2461	}2462 2463	hull = isl_basic_map_finalize(hull);2464	if (input)2465		input->cached_simple_hull[shift] = isl_basic_map_copy(hull);2466	isl_map_free(input);2467 2468	return hull;2469}2470 2471/* Compute a superset of the convex hull of map that is described2472 * by only translates of the constraints in the constituents of map.2473 */2474__isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)2475{2476	return map_simple_hull(map, 1);2477}2478 2479__isl_give isl_basic_set *isl_set_simple_hull(__isl_take isl_set *set)2480{2481	return bset_from_bmap(isl_map_simple_hull(set_to_map(set)));2482}2483 2484/* Compute a superset of the convex hull of map that is described2485 * by only the constraints in the constituents of map.2486 */2487__isl_give isl_basic_map *isl_map_unshifted_simple_hull(2488	__isl_take isl_map *map)2489{2490	return map_simple_hull(map, 0);2491}2492 2493__isl_give isl_basic_set *isl_set_unshifted_simple_hull(2494	__isl_take isl_set *set)2495{2496	return isl_map_unshifted_simple_hull(set);2497}2498 2499/* Drop all inequalities from "bmap1" that do not also appear in "bmap2".2500 * A constraint that appears with different constant terms2501 * in "bmap1" and "bmap2" is also kept, with the least restrictive2502 * (i.e., greatest) constant term.2503 * "bmap1" and "bmap2" are assumed to have the same (known)2504 * integer divisions.2505 * The constraints of both "bmap1" and "bmap2" are assumed2506 * to have been sorted using isl_basic_map_sort_constraints.2507 *2508 * Run through the inequality constraints of "bmap1" and "bmap2"2509 * in sorted order.2510 * Each constraint of "bmap1" without a matching constraint in "bmap2"2511 * is removed.2512 * If a match is found, the constraint is kept.  If needed, the constant2513 * term of the constraint is adjusted.2514 */2515static __isl_give isl_basic_map *select_shared_inequalities(2516	__isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)2517{2518	int i1, i2;2519 2520	bmap1 = isl_basic_map_cow(bmap1);2521	if (!bmap1 || !bmap2)2522		return isl_basic_map_free(bmap1);2523 2524	i1 = bmap1->n_ineq - 1;2525	i2 = bmap2->n_ineq - 1;2526	while (bmap1 && i1 >= 0 && i2 >= 0) {2527		int cmp;2528 2529		cmp = isl_basic_map_constraint_cmp(bmap1, bmap1->ineq[i1],2530							bmap2->ineq[i2]);2531		if (cmp < 0) {2532			--i2;2533			continue;2534		}2535		if (cmp > 0) {2536			if (isl_basic_map_drop_inequality(bmap1, i1) < 0)2537				bmap1 = isl_basic_map_free(bmap1);2538			--i1;2539			continue;2540		}2541		if (isl_int_lt(bmap1->ineq[i1][0], bmap2->ineq[i2][0]))2542			isl_int_set(bmap1->ineq[i1][0], bmap2->ineq[i2][0]);2543		--i1;2544		--i2;2545	}2546	for (; i1 >= 0; --i1)2547		if (isl_basic_map_drop_inequality(bmap1, i1) < 0)2548			bmap1 = isl_basic_map_free(bmap1);2549 2550	return bmap1;2551}2552 2553/* Drop all equalities from "bmap1" that do not also appear in "bmap2".2554 * "bmap1" and "bmap2" are assumed to have the same (known)2555 * integer divisions.2556 *2557 * Run through the equality constraints of "bmap1" and "bmap2".2558 * Each constraint of "bmap1" without a matching constraint in "bmap2"2559 * is removed.2560 */2561static __isl_give isl_basic_map *select_shared_equalities(2562	__isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)2563{2564	int i1, i2;2565	isl_size total;2566 2567	bmap1 = isl_basic_map_cow(bmap1);2568	total = isl_basic_map_dim(bmap1, isl_dim_all);2569	if (total < 0 || !bmap2)2570		return isl_basic_map_free(bmap1);2571 2572	i1 = bmap1->n_eq - 1;2573	i2 = bmap2->n_eq - 1;2574	while (bmap1 && i1 >= 0 && i2 >= 0) {2575		int last1, last2;2576 2577		last1 = isl_seq_last_non_zero(bmap1->eq[i1] + 1, total);2578		last2 = isl_seq_last_non_zero(bmap2->eq[i2] + 1, total);2579		if (last1 > last2) {2580			--i2;2581			continue;2582		}2583		if (last1 < last2) {2584			if (isl_basic_map_drop_equality(bmap1, i1) < 0)2585				bmap1 = isl_basic_map_free(bmap1);2586			--i1;2587			continue;2588		}2589		if (!isl_seq_eq(bmap1->eq[i1], bmap2->eq[i2], 1 + total)) {2590			if (isl_basic_map_drop_equality(bmap1, i1) < 0)2591				bmap1 = isl_basic_map_free(bmap1);2592		}2593		--i1;2594		--i2;2595	}2596	for (; i1 >= 0; --i1)2597		if (isl_basic_map_drop_equality(bmap1, i1) < 0)2598			bmap1 = isl_basic_map_free(bmap1);2599 2600	return bmap1;2601}2602 2603/* Compute a superset of "bmap1" and "bmap2" that is described2604 * by only the constraints that appear in both "bmap1" and "bmap2".2605 *2606 * First drop constraints that involve unknown integer divisions2607 * since it is not trivial to check whether two such integer divisions2608 * in different basic maps are the same.2609 * Then align the remaining (known) divs and sort the constraints.2610 * Finally drop all inequalities and equalities from "bmap1" that2611 * do not also appear in "bmap2".2612 */2613__isl_give isl_basic_map *isl_basic_map_plain_unshifted_simple_hull(2614	__isl_take isl_basic_map *bmap1, __isl_take isl_basic_map *bmap2)2615{2616	if (isl_basic_map_check_equal_space(bmap1, bmap2) < 0)2617		goto error;2618 2619	bmap1 = isl_basic_map_drop_constraints_involving_unknown_divs(bmap1);2620	bmap2 = isl_basic_map_drop_constraints_involving_unknown_divs(bmap2);2621	bmap1 = isl_basic_map_order_divs(bmap1);2622	bmap2 = isl_basic_map_align_divs(bmap2, bmap1);2623	bmap1 = isl_basic_map_align_divs(bmap1, bmap2);2624	bmap1 = isl_basic_map_sort_constraints(bmap1);2625	bmap2 = isl_basic_map_sort_constraints(bmap2);2626 2627	bmap1 = select_shared_inequalities(bmap1, bmap2);2628	bmap1 = select_shared_equalities(bmap1, bmap2);2629 2630	isl_basic_map_free(bmap2);2631	bmap1 = isl_basic_map_finalize(bmap1);2632	return bmap1;2633error:2634	isl_basic_map_free(bmap1);2635	isl_basic_map_free(bmap2);2636	return NULL;2637}2638 2639/* Compute a superset of the convex hull of "map" that is described2640 * by only the constraints in the constituents of "map".2641 * In particular, the result is composed of constraints that appear2642 * in each of the basic maps of "map"2643 *2644 * Constraints that involve unknown integer divisions are dropped2645 * since it is not trivial to check whether two such integer divisions2646 * in different basic maps are the same.2647 *2648 * The hull is initialized from the first basic map and then2649 * updated with respect to the other basic maps in turn.2650 */2651__isl_give isl_basic_map *isl_map_plain_unshifted_simple_hull(2652	__isl_take isl_map *map)2653{2654	int i;2655	isl_basic_map *hull;2656 2657	if (!map)2658		return NULL;2659	if (map->n <= 1)2660		return map_simple_hull_trivial(map);2661	map = isl_map_drop_constraints_involving_unknown_divs(map);2662	hull = isl_basic_map_copy(map->p[0]);2663	for (i = 1; i < map->n; ++i) {2664		isl_basic_map *bmap_i;2665 2666		bmap_i = isl_basic_map_copy(map->p[i]);2667		hull = isl_basic_map_plain_unshifted_simple_hull(hull, bmap_i);2668	}2669 2670	isl_map_free(map);2671	return hull;2672}2673 2674/* Compute a superset of the convex hull of "set" that is described2675 * by only the constraints in the constituents of "set".2676 * In particular, the result is composed of constraints that appear2677 * in each of the basic sets of "set"2678 */2679__isl_give isl_basic_set *isl_set_plain_unshifted_simple_hull(2680	__isl_take isl_set *set)2681{2682	return isl_map_plain_unshifted_simple_hull(set);2683}2684 2685/* Check if "ineq" is a bound on "set" and, if so, add it to "hull".2686 *2687 * For each basic set in "set", we first check if the basic set2688 * contains a translate of "ineq".  If this translate is more relaxed,2689 * then we assume that "ineq" is not a bound on this basic set.2690 * Otherwise, we know that it is a bound.2691 * If the basic set does not contain a translate of "ineq", then2692 * we call is_bound to perform the test.2693 */2694static __isl_give isl_basic_set *add_bound_from_constraint(2695	__isl_take isl_basic_set *hull, struct sh_data *data,2696	__isl_keep isl_set *set, isl_int *ineq)2697{2698	int i, k;2699	isl_ctx *ctx;2700	uint32_t c_hash;2701	struct ineq_cmp_data v;2702	isl_size total;2703 2704	total = isl_basic_set_dim(hull, isl_dim_all);2705	if (total < 0 || !set)2706		return isl_basic_set_free(hull);2707 2708	v.len = total;2709	v.p = ineq;2710	c_hash = isl_seq_get_hash(ineq + 1, v.len);2711 2712	ctx = isl_basic_set_get_ctx(hull);2713	for (i = 0; i < set->n; ++i) {2714		int bound;2715		struct isl_hash_table_entry *entry;2716 2717		entry = isl_hash_table_find(ctx, data->p[i].table,2718						c_hash, &has_ineq, &v, 0);2719		if (!entry)2720			return isl_basic_set_free(hull);2721		if (entry != isl_hash_table_entry_none) {2722			isl_int *ineq_i = entry->data;2723			int neg, more_relaxed;2724 2725			neg = isl_seq_is_neg(ineq_i + 1, ineq + 1, v.len);2726			if (neg)2727				isl_int_neg(ineq_i[0], ineq_i[0]);2728			more_relaxed = isl_int_gt(ineq_i[0], ineq[0]);2729			if (neg)2730				isl_int_neg(ineq_i[0], ineq_i[0]);2731			if (more_relaxed)2732				break;2733			else2734				continue;2735		}2736		bound = is_bound(data, set, i, ineq, 0);2737		if (bound < 0)2738			return isl_basic_set_free(hull);2739		if (!bound)2740			break;2741	}2742	if (i < set->n)2743		return hull;2744 2745	k = isl_basic_set_alloc_inequality(hull);2746	if (k < 0)2747		return isl_basic_set_free(hull);2748	isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);2749 2750	return hull;2751}2752 2753/* Compute a superset of the convex hull of "set" that is described2754 * by only some of the "n_ineq" constraints in the list "ineq", where "set"2755 * has no parameters or integer divisions.2756 *2757 * The inequalities in "ineq" are assumed to have been sorted such2758 * that constraints with the same linear part appear together and2759 * that among constraints with the same linear part, those with2760 * smaller constant term appear first.2761 *2762 * We reuse the same data structure that is used by uset_simple_hull,2763 * but we do not need the hull table since we will not consider the2764 * same constraint more than once.  We therefore allocate it with zero size.2765 *2766 * We run through the constraints and try to add them one by one,2767 * skipping identical constraints.  If we have added a constraint and2768 * the next constraint is a more relaxed translate, then we skip this2769 * next constraint as well.2770 */2771static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_constraints(2772	__isl_take isl_set *set, int n_ineq, isl_int **ineq)2773{2774	int i;2775	int last_added = 0;2776	struct sh_data *data = NULL;2777	isl_basic_set *hull = NULL;2778	isl_size dim;2779 2780	hull = isl_basic_set_alloc_space(isl_set_get_space(set), 0, 0, n_ineq);2781	if (!hull)2782		goto error;2783 2784	data = sh_data_alloc(set, 0);2785	if (!data)2786		goto error;2787 2788	dim = isl_set_dim(set, isl_dim_set);2789	if (dim < 0)2790		goto error;2791	for (i = 0; i < n_ineq; ++i) {2792		int hull_n_ineq = hull->n_ineq;2793		int parallel;2794 2795		parallel = i > 0 && isl_seq_eq(ineq[i - 1] + 1, ineq[i] + 1,2796						dim);2797		if (parallel &&2798		    (last_added || isl_int_eq(ineq[i - 1][0], ineq[i][0])))2799			continue;2800		hull = add_bound_from_constraint(hull, data, set, ineq[i]);2801		if (!hull)2802			goto error;2803		last_added = hull->n_ineq > hull_n_ineq;2804	}2805 2806	sh_data_free(data);2807	isl_set_free(set);2808	return hull;2809error:2810	sh_data_free(data);2811	isl_set_free(set);2812	isl_basic_set_free(hull);2813	return NULL;2814}2815 2816/* Collect pointers to all the inequalities in the elements of "list"2817 * in "ineq".  For equalities, store both a pointer to the equality and2818 * a pointer to its opposite, which is first copied to "mat".2819 * "ineq" and "mat" are assumed to have been preallocated to the right size2820 * (the number of inequalities + 2 times the number of equalites and2821 * the number of equalities, respectively).2822 */2823static __isl_give isl_mat *collect_inequalities(__isl_take isl_mat *mat,2824	__isl_keep isl_basic_set_list *list, isl_int **ineq)2825{2826	int i, j, n_eq, n_ineq;2827	isl_size n;2828 2829	n = isl_basic_set_list_n_basic_set(list);2830	if (!mat || n < 0)2831		return isl_mat_free(mat);2832 2833	n_eq = 0;2834	n_ineq = 0;2835	for (i = 0; i < n; ++i) {2836		isl_basic_set *bset;2837		bset = isl_basic_set_list_get_basic_set(list, i);2838		if (!bset)2839			return isl_mat_free(mat);2840		for (j = 0; j < bset->n_eq; ++j) {2841			ineq[n_ineq++] = mat->row[n_eq];2842			ineq[n_ineq++] = bset->eq[j];2843			isl_seq_neg(mat->row[n_eq++], bset->eq[j], mat->n_col);2844		}2845		for (j = 0; j < bset->n_ineq; ++j)2846			ineq[n_ineq++] = bset->ineq[j];2847		isl_basic_set_free(bset);2848	}2849 2850	return mat;2851}2852 2853/* Comparison routine for use as an isl_sort callback.2854 *2855 * Constraints with the same linear part are sorted together and2856 * among constraints with the same linear part, those with smaller2857 * constant term are sorted first.2858 */2859static int cmp_ineq(const void *a, const void *b, void *arg)2860{2861	unsigned dim = *(unsigned *) arg;2862	isl_int * const *ineq1 = a;2863	isl_int * const *ineq2 = b;2864	int cmp;2865 2866	cmp = isl_seq_cmp((*ineq1) + 1, (*ineq2) + 1, dim);2867	if (cmp != 0)2868		return cmp;2869	return isl_int_cmp((*ineq1)[0], (*ineq2)[0]);2870}2871 2872/* Compute a superset of the convex hull of "set" that is described2873 * by only constraints in the elements of "list", where "set" has2874 * no parameters or integer divisions.2875 *2876 * We collect all the constraints in those elements and then2877 * sort the constraints such that constraints with the same linear part2878 * are sorted together and that those with smaller constant term are2879 * sorted first.2880 */2881static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_basic_set_list(2882	__isl_take isl_set *set, __isl_take isl_basic_set_list *list)2883{2884	int i, n_eq, n_ineq;2885	isl_size n;2886	isl_size dim;2887	isl_ctx *ctx;2888	isl_mat *mat = NULL;2889	isl_int **ineq = NULL;2890	isl_basic_set *hull;2891 2892	n = isl_basic_set_list_n_basic_set(list);2893	if (!set || n < 0)2894		goto error;2895	ctx = isl_set_get_ctx(set);2896 2897	n_eq = 0;2898	n_ineq = 0;2899	for (i = 0; i < n; ++i) {2900		isl_basic_set *bset;2901		bset = isl_basic_set_list_get_basic_set(list, i);2902		if (!bset)2903			goto error;2904		n_eq += bset->n_eq;2905		n_ineq += 2 * bset->n_eq + bset->n_ineq;2906		isl_basic_set_free(bset);2907	}2908 2909	ineq = isl_alloc_array(ctx, isl_int *, n_ineq);2910	if (n_ineq > 0 && !ineq)2911		goto error;2912 2913	dim = isl_set_dim(set, isl_dim_set);2914	if (dim < 0)2915		goto error;2916	mat = isl_mat_alloc(ctx, n_eq, 1 + dim);2917	mat = collect_inequalities(mat, list, ineq);2918	if (!mat)2919		goto error;2920 2921	if (isl_sort(ineq, n_ineq, sizeof(ineq[0]), &cmp_ineq, &dim) < 0)2922		goto error;2923 2924	hull = uset_unshifted_simple_hull_from_constraints(set, n_ineq, ineq);2925 2926	isl_mat_free(mat);2927	free(ineq);2928	isl_basic_set_list_free(list);2929	return hull;2930error:2931	isl_mat_free(mat);2932	free(ineq);2933	isl_set_free(set);2934	isl_basic_set_list_free(list);2935	return NULL;2936}2937 2938/* Compute a superset of the convex hull of "map" that is described2939 * by only constraints in the elements of "list".2940 *2941 * If the list is empty, then we can only describe the universe set.2942 * If the input map is empty, then all constraints are valid, so2943 * we return the intersection of the elements in "list".2944 *2945 * Otherwise, we align all divs and temporarily treat them2946 * as regular variables, computing the unshifted simple hull in2947 * uset_unshifted_simple_hull_from_basic_set_list.2948 */2949static __isl_give isl_basic_map *map_unshifted_simple_hull_from_basic_map_list(2950	__isl_take isl_map *map, __isl_take isl_basic_map_list *list)2951{2952	isl_size n;2953	isl_basic_map *model;2954	isl_basic_map *hull;2955	isl_set *set;2956	isl_basic_set_list *bset_list;2957 2958	n = isl_basic_map_list_n_basic_map(list);2959	if (!map || n < 0)2960		goto error;2961 2962	if (n == 0) {2963		isl_space *space;2964 2965		space = isl_map_get_space(map);2966		isl_map_free(map);2967		isl_basic_map_list_free(list);2968		return isl_basic_map_universe(space);2969	}2970	if (isl_map_plain_is_empty(map)) {2971		isl_map_free(map);2972		return isl_basic_map_list_intersect(list);2973	}2974 2975	map = isl_map_align_divs_to_basic_map_list(map, list);2976	if (!map)2977		goto error;2978	list = isl_basic_map_list_align_divs_to_basic_map(list, map->p[0]);2979 2980	model = isl_basic_map_list_get_basic_map(list, 0);2981 2982	set = isl_map_underlying_set(map);2983	bset_list = isl_basic_map_list_underlying_set(list);2984 2985	hull = uset_unshifted_simple_hull_from_basic_set_list(set, bset_list);2986	hull = isl_basic_map_overlying_set(hull, model);2987 2988	return hull;2989error:2990	isl_map_free(map);2991	isl_basic_map_list_free(list);2992	return NULL;2993}2994 2995/* Return a sequence of the basic maps that make up the maps in "list".2996 */2997static __isl_give isl_basic_map_list *collect_basic_maps(2998	__isl_take isl_map_list *list)2999{3000	int i;3001	isl_size n;3002	isl_ctx *ctx;3003	isl_basic_map_list *bmap_list;3004 3005	if (!list)3006		return NULL;3007	n = isl_map_list_n_map(list);3008	ctx = isl_map_list_get_ctx(list);3009	bmap_list = isl_basic_map_list_alloc(ctx, 0);3010	if (n < 0)3011		bmap_list = isl_basic_map_list_free(bmap_list);3012 3013	for (i = 0; i < n; ++i) {3014		isl_map *map;3015		isl_basic_map_list *list_i;3016 3017		map = isl_map_list_get_map(list, i);3018		map = isl_map_compute_divs(map);3019		list_i = isl_map_get_basic_map_list(map);3020		isl_map_free(map);3021		bmap_list = isl_basic_map_list_concat(bmap_list, list_i);3022	}3023 3024	isl_map_list_free(list);3025	return bmap_list;3026}3027 3028/* Compute a superset of the convex hull of "map" that is described3029 * by only constraints in the elements of "list".3030 *3031 * If "map" is the universe, then the convex hull (and therefore3032 * any superset of the convexhull) is the universe as well.3033 *3034 * Otherwise, we collect all the basic maps in the map list and3035 * continue with map_unshifted_simple_hull_from_basic_map_list.3036 */3037__isl_give isl_basic_map *isl_map_unshifted_simple_hull_from_map_list(3038	__isl_take isl_map *map, __isl_take isl_map_list *list)3039{3040	isl_basic_map_list *bmap_list;3041	int is_universe;3042 3043	is_universe = isl_map_plain_is_universe(map);3044	if (is_universe < 0)3045		map = isl_map_free(map);3046	if (is_universe < 0 || is_universe) {3047		isl_map_list_free(list);3048		return isl_map_unshifted_simple_hull(map);3049	}3050 3051	bmap_list = collect_basic_maps(list);3052	return map_unshifted_simple_hull_from_basic_map_list(map, bmap_list);3053}3054 3055/* Compute a superset of the convex hull of "set" that is described3056 * by only constraints in the elements of "list".3057 */3058__isl_give isl_basic_set *isl_set_unshifted_simple_hull_from_set_list(3059	__isl_take isl_set *set, __isl_take isl_set_list *list)3060{3061	return isl_map_unshifted_simple_hull_from_map_list(set, list);3062}3063 3064/* Given a set "set", return parametric bounds on the dimension "dim".3065 */3066static __isl_give isl_basic_set *set_bounds(__isl_keep isl_set *set, int dim)3067{3068	isl_size set_dim = isl_set_dim(set, isl_dim_set);3069	if (set_dim < 0)3070		return NULL;3071	set = isl_set_copy(set);3072	set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1));3073	set = isl_set_eliminate_dims(set, 0, dim);3074	return isl_set_convex_hull(set);3075}3076 3077/* Computes a "simple hull" and then check if each dimension in the3078 * resulting hull is bounded by a symbolic constant.  If not, the3079 * hull is intersected with the corresponding bounds on the whole set.3080 */3081__isl_give isl_basic_set *isl_set_bounded_simple_hull(__isl_take isl_set *set)3082{3083	int i, j;3084	struct isl_basic_set *hull;3085	isl_size nparam, dim, total;3086	unsigned left;3087	int removed_divs = 0;3088 3089	hull = isl_set_simple_hull(isl_set_copy(set));3090	nparam = isl_basic_set_dim(hull, isl_dim_param);3091	dim = isl_basic_set_dim(hull, isl_dim_set);3092	total = isl_basic_set_dim(hull, isl_dim_all);3093	if (nparam < 0 || dim < 0 || total < 0)3094		goto error;3095 3096	for (i = 0; i < dim; ++i) {3097		int lower = 0, upper = 0;3098		struct isl_basic_set *bounds;3099 3100		left = total - nparam - i - 1;3101		for (j = 0; j < hull->n_eq; ++j) {3102			if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))3103				continue;3104			if (isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1,3105						    left) == -1)3106				break;3107		}3108		if (j < hull->n_eq)3109			continue;3110 3111		for (j = 0; j < hull->n_ineq; ++j) {3112			if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))3113				continue;3114			if (isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1,3115						    left) != -1 ||3116			    isl_seq_first_non_zero(hull->ineq[j]+1+nparam,3117						    i) != -1)3118				continue;3119			if (isl_int_is_pos(hull->ineq[j][1 + nparam + i]))3120				lower = 1;3121			else3122				upper = 1;3123			if (lower && upper)3124				break;3125		}3126 3127		if (lower && upper)3128			continue;3129 3130		if (!removed_divs) {3131			set = isl_set_remove_divs(set);3132			if (!set)3133				goto error;3134			removed_divs = 1;3135		}3136		bounds = set_bounds(set, i);3137		hull = isl_basic_set_intersect(hull, bounds);3138		if (!hull)3139			goto error;3140	}3141 3142	isl_set_free(set);3143	return hull;3144error:3145	isl_set_free(set);3146	isl_basic_set_free(hull);3147	return NULL;3148}3149