/*

 * Copyright 2008-2009 Katholieke Universiteit Leuven

 * Copyright 2010      INRIA Saclay

 * Copyright 2012-2013 Ecole Normale Superieure

 *

 * Use of this software is governed by the MIT license

 *

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

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

 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,

 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France 

 * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France

 */

#include "isl_map_private.h"

#include <isl_seq.h>

#include <isl/options.h>

#include "isl_tab.h"

#include <isl_mat_private.h>

#include <isl_local_space_private.h>

#include <isl_vec_private.h>

#define STATUS_ERROR		-1

#define STATUS_REDUNDANT	 1

#define STATUS_VALID	 	 2

#define STATUS_SEPARATE	 	 3

#define STATUS_CUT	 	 4

#define STATUS_ADJ_EQ	 	 5

#define STATUS_ADJ_INEQ	 	 6

static int status_in(isl_int *ineq, struct isl_tab *tab)

{

	enum isl_ineq_type type = isl_tab_ineq_type(tab, ineq);

	switch (type) {

	default:

	case isl_ineq_error:		return STATUS_ERROR;

	case isl_ineq_redundant:	return STATUS_VALID;

	case isl_ineq_separate:		return STATUS_SEPARATE;

	case isl_ineq_cut:		return STATUS_CUT;

	case isl_ineq_adj_eq:		return STATUS_ADJ_EQ;

	case isl_ineq_adj_ineq:		return STATUS_ADJ_INEQ;

	}

}

/* Compute the position of the equalities of basic map "bmap_i"

 * with respect to the basic map represented by "tab_j".

 * The resulting array has twice as many entries as the number

 * of equalities corresponding to the two inequalties to which

 * each equality corresponds.

 */

static int *eq_status_in(__isl_keep isl_basic_map *bmap_i,

	struct isl_tab *tab_j)

{

	int k, l;

	int *eq = isl_calloc_array(bmap_i->ctx, int, 2 * bmap_i->n_eq);

	unsigned dim;

	if (!eq)

		return NULL;

	dim = isl_basic_map_total_dim(bmap_i);

	for (k = 0; k < bmap_i->n_eq; ++k) {

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

			isl_seq_neg(bmap_i->eq[k], bmap_i->eq[k], 1+dim);

			eq[2 * k + l] = status_in(bmap_i->eq[k], tab_j);

			if (eq[2 * k + l] == STATUS_ERROR)

				goto error;

		}

		if (eq[2 * k] == STATUS_SEPARATE ||

		    eq[2 * k + 1] == STATUS_SEPARATE)

			break;

	}

	return eq;

error:

	free(eq);

	return NULL;

}

/* Compute the position of the inequalities of basic map "bmap_i"

 * (also represented by "tab_i", if not NULL) with respect to the basic map

 * represented by "tab_j".

 */

static int *ineq_status_in(__isl_keep isl_basic_map *bmap_i,

	struct isl_tab *tab_i, struct isl_tab *tab_j)

{

	int k;

	unsigned n_eq = bmap_i->n_eq;

	int *ineq = isl_calloc_array(bmap_i->ctx, int, bmap_i->n_ineq);

	if (!ineq)

		return NULL;

	for (k = 0; k < bmap_i->n_ineq; ++k) {

		if (tab_i && isl_tab_is_redundant(tab_i, n_eq + k)) {

			ineq[k] = STATUS_REDUNDANT;

			continue;

		}

		ineq[k] = status_in(bmap_i->ineq[k], tab_j);

		if (ineq[k] == STATUS_ERROR)

			goto error;

		if (ineq[k] == STATUS_SEPARATE)

			break;

	}

	return ineq;

error:

	free(ineq);

	return NULL;

}

static int any(int *con, unsigned len, int status)

{

	int i;

	for (i = 0; i < len ; ++i)

		if (con[i] == status)

			return 1;

	return 0;

}

static int count(int *con, unsigned len, int status)

{

	int i;

	int c = 0;

	for (i = 0; i < len ; ++i)

		if (con[i] == status)

			c++;

	return c;

}

static int all(int *con, unsigned len, int status)

{

	int i;

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

		if (con[i] == STATUS_REDUNDANT)

			continue;

		if (con[i] != status)

			return 0;

	}

	return 1;

}

static void drop(struct isl_map *map, int i, struct isl_tab **tabs)

{

	isl_basic_map_free(map->p[i]);

	isl_tab_free(tabs[i]);

	if (i != map->n - 1) {

		map->p[i] = map->p[map->n - 1];

		tabs[i] = tabs[map->n - 1];

	}

	tabs[map->n - 1] = NULL;

	map->n--;

}

/* Replace the pair of basic maps i and j by the basic map bounded

 * by the valid constraints in both basic maps and the constraint

 * in extra (if not NULL).

 */

static int fuse(struct isl_map *map, int i, int j,

	struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j,

	__isl_keep isl_mat *extra)

{

	int k, l;

	struct isl_basic_map *fused = NULL;

	struct isl_tab *fused_tab = NULL;

	unsigned total = isl_basic_map_total_dim(map->p[i]);

	unsigned extra_rows = extra ? extra->n_row : 0;

	fused = isl_basic_map_alloc_space(isl_space_copy(map->p[i]->dim),

			map->p[i]->n_div,

			map->p[i]->n_eq + map->p[j]->n_eq,

			map->p[i]->n_ineq + map->p[j]->n_ineq + extra_rows);

	if (!fused)

		goto error;

	for (k = 0; k < map->p[i]->n_eq; ++k) {

		if (eq_i && (eq_i[2 * k] != STATUS_VALID ||

			     eq_i[2 * k + 1] != STATUS_VALID))

			continue;

		l = isl_basic_map_alloc_equality(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->eq[l], map->p[i]->eq[k], 1 + total);

	}

	for (k = 0; k < map->p[j]->n_eq; ++k) {

		if (eq_j && (eq_j[2 * k] != STATUS_VALID ||

			     eq_j[2 * k + 1] != STATUS_VALID))

			continue;

		l = isl_basic_map_alloc_equality(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->eq[l], map->p[j]->eq[k], 1 + total);

	}

	for (k = 0; k < map->p[i]->n_ineq; ++k) {

		if (ineq_i[k] != STATUS_VALID)

			continue;

		l = isl_basic_map_alloc_inequality(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->ineq[l], map->p[i]->ineq[k], 1 + total);

	}

	for (k = 0; k < map->p[j]->n_ineq; ++k) {

		if (ineq_j[k] != STATUS_VALID)

			continue;

		l = isl_basic_map_alloc_inequality(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->ineq[l], map->p[j]->ineq[k], 1 + total);

	}

	for (k = 0; k < map->p[i]->n_div; ++k) {

		int l = isl_basic_map_alloc_div(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->div[l], map->p[i]->div[k], 1 + 1 + total);

	}

	for (k = 0; k < extra_rows; ++k) {

		l = isl_basic_map_alloc_inequality(fused);

		if (l < 0)

			goto error;

		isl_seq_cpy(fused->ineq[l], extra->row[k], 1 + total);

	}

	fused = isl_basic_map_gauss(fused, NULL);

	ISL_F_SET(fused, ISL_BASIC_MAP_FINAL);

	if (ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_RATIONAL) &&

	    ISL_F_ISSET(map->p[j], ISL_BASIC_MAP_RATIONAL))

		ISL_F_SET(fused, ISL_BASIC_MAP_RATIONAL);

	fused_tab = isl_tab_from_basic_map(fused, 0);

	if (isl_tab_detect_redundant(fused_tab) < 0)

		goto error;

	isl_basic_map_free(map->p[i]);

	map->p[i] = fused;

	isl_tab_free(tabs[i]);

	tabs[i] = fused_tab;

	drop(map, j, tabs);

	return 1;

error:

	isl_tab_free(fused_tab);

	isl_basic_map_free(fused);

	return -1;

}

/* Given a pair of basic maps i and j such that all constraints are either

 * "valid" or "cut", check if the facets corresponding to the "cut"

 * constraints of i lie entirely within basic map j.

 * If so, replace the pair by the basic map consisting of the valid

 * constraints in both basic maps.

 *

 * To see that we are not introducing any extra points, call the

 * two basic maps A and B and the resulting map U and let x

 * be an element of U \setminus ( A \cup B ).

 * Then there is a pair of cut constraints c_1 and c_2 in A and B such that x

 * violates them.  Let X be the intersection of U with the opposites

 * of these constraints.  Then x \in X.

 * The facet corresponding to c_1 contains the corresponding facet of A.

 * This facet is entirely contained in B, so c_2 is valid on the facet.

 * However, since it is also (part of) a facet of X, -c_2 is also valid

 * on the facet.  This means c_2 is saturated on the facet, so c_1 and

 * c_2 must be opposites of each other, but then x could not violate

 * both of them.

 */

static int check_facets(struct isl_map *map, int i, int j,

	struct isl_tab **tabs, int *ineq_i, int *ineq_j)

{

	int k, l;

	struct isl_tab_undo *snap;

	unsigned n_eq = map->p[i]->n_eq;

	snap = isl_tab_snap(tabs[i]);

	for (k = 0; k < map->p[i]->n_ineq; ++k) {

		if (ineq_i[k] != STATUS_CUT)

			continue;

		if (isl_tab_select_facet(tabs[i], n_eq + k) < 0)

			return -1;

		for (l = 0; l < map->p[j]->n_ineq; ++l) {

			int stat;

			if (ineq_j[l] != STATUS_CUT)

				continue;

			stat = status_in(map->p[j]->ineq[l], tabs[i]);

			if (stat != STATUS_VALID)

				break;

		}

		if (isl_tab_rollback(tabs[i], snap) < 0)

			return -1;

		if (l < map->p[j]->n_ineq)

			break;

	}

	if (k < map->p[i]->n_ineq)

		/* BAD CUT PAIR */

		return 0;

	return fuse(map, i, j, tabs, NULL, ineq_i, NULL, ineq_j, NULL);

}

/* Check if basic map "i" contains the basic map represented

 * by the tableau "tab".

 */

static int contains(struct isl_map *map, int i, int *ineq_i,

	struct isl_tab *tab)

{

	int k, l;

	unsigned dim;

	dim = isl_basic_map_total_dim(map->p[i]);

	for (k = 0; k < map->p[i]->n_eq; ++k) {

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

			int stat;

			isl_seq_neg(map->p[i]->eq[k], map->p[i]->eq[k], 1+dim);

			stat = status_in(map->p[i]->eq[k], tab);

			if (stat != STATUS_VALID)

				return 0;

		}

	}

	for (k = 0; k < map->p[i]->n_ineq; ++k) {

		int stat;

		if (ineq_i[k] == STATUS_REDUNDANT)

			continue;

		stat = status_in(map->p[i]->ineq[k], tab);

		if (stat != STATUS_VALID)

			return 0;

	}

	return 1;

}

/* Basic map "i" has an inequality (say "k") that is adjacent

 * to some inequality of basic map "j".  All the other inequalities

 * are valid for "j".

 * Check if basic map "j" forms an extension of basic map "i".

 *

 * Note that this function is only called if some of the equalities or

 * inequalities of basic map "j" do cut basic map "i".  The function is

 * correct even if there are no such cut constraints, but in that case

 * the additional checks performed by this function are overkill.

 *

 * In particular, we replace constraint k, say f >= 0, by constraint

 * f <= -1, add the inequalities of "j" that are valid for "i"

 * and check if the result is a subset of basic map "j".

 * If so, then we know that this result is exactly equal to basic map "j"

 * since all its constraints are valid for basic map "j".

 * By combining the valid constraints of "i" (all equalities and all

 * inequalities except "k") and the valid constraints of "j" we therefore

 * obtain a basic map that is equal to their union.

 * In this case, there is no need to perform a rollback of the tableau

 * since it is going to be destroyed in fuse().

 *

 *

 *	|\__			|\__

 *	|   \__			|   \__

 *	|      \_	=>	|      \__

 *	|_______| _		|_________\

 *

 *

 *	|\			|\

 *	| \			| \

 *	|  \			|  \

 *	|  |			|   \

 *	|  ||\		=>      |    \

 *	|  || \			|     \

 *	|  ||  |		|      |

 *	|__||_/			|_____/

 */

static int is_adj_ineq_extension(__isl_keep isl_map *map, int i, int j,

	struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)

{

	int k;

	struct isl_tab_undo *snap;

	unsigned n_eq = map->p[i]->n_eq;

	unsigned total = isl_basic_map_total_dim(map->p[i]);

	int r;

	if (isl_tab_extend_cons(tabs[i], 1 + map->p[j]->n_ineq) < 0)

		return -1;

	for (k = 0; k < map->p[i]->n_ineq; ++k)

		if (ineq_i[k] == STATUS_ADJ_INEQ)

			break;

	if (k >= map->p[i]->n_ineq)

		isl_die(isl_map_get_ctx(map), isl_error_internal,

			"ineq_i should have exactly one STATUS_ADJ_INEQ",

			return -1);

	snap = isl_tab_snap(tabs[i]);

	if (isl_tab_unrestrict(tabs[i], n_eq + k) < 0)

		return -1;

	isl_seq_neg(map->p[i]->ineq[k], map->p[i]->ineq[k], 1 + total);

	isl_int_sub_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);

	r = isl_tab_add_ineq(tabs[i], map->p[i]->ineq[k]);

	isl_seq_neg(map->p[i]->ineq[k], map->p[i]->ineq[k], 1 + total);

	isl_int_sub_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);

	if (r < 0)

		return -1;

	for (k = 0; k < map->p[j]->n_ineq; ++k) {

		if (ineq_j[k] != STATUS_VALID)

			continue;

		if (isl_tab_add_ineq(tabs[i], map->p[j]->ineq[k]) < 0)

			return -1;

	}

	if (contains(map, j, ineq_j, tabs[i]))

		return fuse(map, i, j, tabs, eq_i, ineq_i, eq_j, ineq_j, NULL);

	if (isl_tab_rollback(tabs[i], snap) < 0)

		return -1;

	return 0;

}

/* Both basic maps have at least one inequality with and adjacent

 * (but opposite) inequality in the other basic map.

 * Check that there are no cut constraints and that there is only

 * a single pair of adjacent inequalities.

 * If so, we can replace the pair by a single basic map described

 * by all but the pair of adjacent inequalities.

 * Any additional points introduced lie strictly between the two

 * adjacent hyperplanes and can therefore be integral.

 *

 *        ____			  _____

 *       /    ||\		 /     \

 *      /     || \		/       \

 *      \     ||  \	=>	\        \

 *       \    ||  /		 \       /

 *        \___||_/		  \_____/

 *

 * The test for a single pair of adjancent inequalities is important

 * for avoiding the combination of two basic maps like the following

 *

 *       /|

 *      / |

 *     /__|

 *         _____

 *         |   |

 *         |   |

 *         |___|

 *

 * If there are some cut constraints on one side, then we may

 * still be able to fuse the two basic maps, but we need to perform

 * some additional checks in is_adj_ineq_extension.

 */

static int check_adj_ineq(struct isl_map *map, int i, int j,

	struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)

{

	int count_i, count_j;

	int cut_i, cut_j;

	count_i = count(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ);

	count_j = count(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ);

	if (count_i != 1 && count_j != 1)

		return 0;

	cut_i = any(eq_i, 2 * map->p[i]->n_eq, STATUS_CUT) ||

		any(ineq_i, map->p[i]->n_ineq, STATUS_CUT);

	cut_j = any(eq_j, 2 * map->p[j]->n_eq, STATUS_CUT) ||

		any(ineq_j, map->p[j]->n_ineq, STATUS_CUT);

	if (!cut_i && !cut_j && count_i == 1 && count_j == 1)

		return fuse(map, i, j, tabs, NULL, ineq_i, NULL, ineq_j, NULL);

	if (count_i == 1 && !cut_i)

		return is_adj_ineq_extension(map, i, j, tabs,

						eq_i, ineq_i, eq_j, ineq_j);

	if (count_j == 1 && !cut_j)

		return is_adj_ineq_extension(map, j, i, tabs,

						eq_j, ineq_j, eq_i, ineq_i);

	return 0;

}

/* Basic map "i" has an inequality "k" that is adjacent to some equality

 * of basic map "j".  All the other inequalities are valid for "j".

 * Check if basic map "j" forms an extension of basic map "i".

 *

 * In particular, we relax constraint "k", compute the corresponding

 * facet and check whether it is included in the other basic map.

 * If so, we know that relaxing the constraint extends the basic

 * map with exactly the other basic map (we already know that this

 * other basic map is included in the extension, because there

 * were no "cut" inequalities in "i") and we can replace the

 * two basic maps by this extension.

 *        ____			  _____

 *       /    || 		 /     |

 *      /     ||  		/      |

 *      \     ||   	=>	\      |

 *       \    ||		 \     |

 *        \___||		  \____|

 */

static int is_adj_eq_extension(struct isl_map *map, int i, int j, int k,

	struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)

{

	int changed = 0;

	int super;

	struct isl_tab_undo *snap, *snap2;

	unsigned n_eq = map->p[i]->n_eq;

	if (isl_tab_is_equality(tabs[i], n_eq + k))

		return 0;

	snap = isl_tab_snap(tabs[i]);

	tabs[i] = isl_tab_relax(tabs[i], n_eq + k);

	snap2 = isl_tab_snap(tabs[i]);

	if (isl_tab_select_facet(tabs[i], n_eq + k) < 0)

		return -1;

	super = contains(map, j, ineq_j, tabs[i]);

	if (super) {

		if (isl_tab_rollback(tabs[i], snap2) < 0)

			return -1;

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

		if (!map->p[i])

			return -1;

		isl_int_add_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);

		ISL_F_SET(map->p[i], ISL_BASIC_MAP_FINAL);

		drop(map, j, tabs);

		changed = 1;

	} else

		if (isl_tab_rollback(tabs[i], snap) < 0)

			return -1;

	return changed;

}

/* Data structure that keeps track of the wrapping constraints

 * and of information to bound the coefficients of those constraints.

 *

 * bound is set if we want to apply a bound on the coefficients

 * mat contains the wrapping constraints

 * max is the bound on the coefficients (if bound is set)

 */

struct isl_wraps {

	int bound;

	isl_mat *mat;

	isl_int max;

};

/* Update wraps->max to be greater than or equal to the coefficients

 * in the equalities and inequalities of bmap that can be removed if we end up

 * applying wrapping.

 */

static void wraps_update_max(struct isl_wraps *wraps,

	__isl_keep isl_basic_map *bmap, int *eq, int *ineq)

{

	int k;

	isl_int max_k;

	unsigned total = isl_basic_map_total_dim(bmap);

	isl_int_init(max_k);

	for (k = 0; k < bmap->n_eq; ++k) {

		if (eq[2 * k] == STATUS_VALID &&

		    eq[2 * k + 1] == STATUS_VALID)

			continue;

		isl_seq_abs_max(bmap->eq[k] + 1, total, &max_k);

		if (isl_int_abs_gt(max_k, wraps->max))

			isl_int_set(wraps->max, max_k);

	}

	for (k = 0; k < bmap->n_ineq; ++k) {

		if (ineq[k] == STATUS_VALID || ineq[k] == STATUS_REDUNDANT)

			continue;

		isl_seq_abs_max(bmap->ineq[k] + 1, total, &max_k);

		if (isl_int_abs_gt(max_k, wraps->max))

			isl_int_set(wraps->max, max_k);

	}

	isl_int_clear(max_k);

}

/* Initialize the isl_wraps data structure.

 * If we want to bound the coefficients of the wrapping constraints,

 * we set wraps->max to the largest coefficient

 * in the equalities and inequalities that can be removed if we end up

 * applying wrapping.

 */

static void wraps_init(struct isl_wraps *wraps, __isl_take isl_mat *mat,

	__isl_keep isl_map *map, int i, int j,

	int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)

{

	isl_ctx *ctx;

	wraps->bound = 0;

	wraps->mat = mat;

	if (!mat)

		return;

	ctx = isl_mat_get_ctx(mat);

	wraps->bound = isl_options_get_coalesce_bounded_wrapping(ctx);

	if (!wraps->bound)

		return;

	isl_int_init(wraps->max);

	isl_int_set_si(wraps->max, 0);

	wraps_update_max(wraps, map->p[i], eq_i, ineq_i);

	wraps_update_max(wraps, map->p[j], eq_j, ineq_j);

}

/* Free the contents of the isl_wraps data structure.

 */

static void wraps_free(struct isl_wraps *wraps)

{

	isl_mat_free(wraps->mat);

	if (wraps->bound)

		isl_int_clear(wraps->max);

}

/* Is the wrapping constraint in row "row" allowed?

 *

 * If wraps->bound is set, we check that none of the coefficients

 * is greater than wraps->max.

 */

static int allow_wrap(struct isl_wraps *wraps, int row)

{

	int i;

	if (!wraps->bound)

		return 1;

	for (i = 1; i < wraps->mat->n_col; ++i)

		if (isl_int_abs_gt(wraps->mat->row[row][i], wraps->max))

			return 0;

	return 1;

}

/* For each non-redundant constraint in "bmap" (as determined by "tab"),

 * wrap the constraint around "bound" such that it includes the whole

 * set "set" and append the resulting constraint to "wraps".

 * "wraps" is assumed to have been pre-allocated to the appropriate size.

 * wraps->n_row is the number of actual wrapped constraints that have

 * been added.

 * If any of the wrapping problems results in a constraint that is

 * identical to "bound", then this means that "set" is unbounded in such

 * way that no wrapping is possible.  If this happens then wraps->n_row

 * is reset to zero.

 * Similarly, if we want to bound the coefficients of the wrapping

 * constraints and a newly added wrapping constraint does not

 * satisfy the bound, then wraps->n_row is also reset to zero.

 */

static int add_wraps(struct isl_wraps *wraps, __isl_keep isl_basic_map *bmap,

	struct isl_tab *tab, isl_int *bound, __isl_keep isl_set *set)

{

	int l;

	int w;

	unsigned total = isl_basic_map_total_dim(bmap);

	w = wraps->mat->n_row;

	for (l = 0; l < bmap->n_ineq; ++l) {

		if (isl_seq_is_neg(bound, bmap->ineq[l], 1 + total))

			continue;

		if (isl_seq_eq(bound, bmap->ineq[l], 1 + total))

			continue;

		if (isl_tab_is_redundant(tab, bmap->n_eq + l))

			continue;

		isl_seq_cpy(wraps->mat->row[w], bound, 1 + total);

		if (!isl_set_wrap_facet(set, wraps->mat->row[w], bmap->ineq[l]))

			return -1;

		if (isl_seq_eq(wraps->mat->row[w], bound, 1 + total))

			goto unbounded;

		if (!allow_wrap(wraps, w))

			goto unbounded;

		++w;

	}

	for (l = 0; l < bmap->n_eq; ++l) {

		if (isl_seq_is_neg(bound, bmap->eq[l], 1 + total))

			continue;

		if (isl_seq_eq(bound, bmap->eq[l], 1 + total))

			continue;

		isl_seq_cpy(wraps->mat->row[w], bound, 1 + total);

		isl_seq_neg(wraps->mat->row[w + 1], bmap->eq[l], 1 + total);

		if (!isl_set_wrap_facet(set, wraps->mat->row[w],

					wraps->mat->row[w + 1]))

			return -1;

		if (isl_seq_eq(wraps->mat->row[w], bound, 1 + total))

			goto unbounded;

		if (!allow_wrap(wraps, w))

			goto unbounded;

		++w;

		isl_seq_cpy(wraps->mat->row[w], bound, 1 + total);

		if (!isl_set_wrap_facet(set, wraps->mat->row[w], bmap->eq[l]))

			return -1;

		if (isl_seq_eq(wraps->mat->row[w], bound, 1 + total))

			goto unbounded;

		if (!allow_wrap(wraps, w))

			goto unbounded;

		++w;

	}

	wraps->mat->n_row = w;

	return 0;

unbounded:

	wraps->mat->n_row = 0;

	return 0;

}

/* Check if the constraints in "wraps" from "first" until the last

 * are all valid for the basic set represented by "tab".

 * If not, wraps->n_row is set to zero.

 */

static int check_wraps(__isl_keep isl_mat *wraps, int first,

	struct isl_tab *tab)

{

	int i;

	for (i = first; i < wraps->n_row; ++i) {

		enum isl_ineq_type type;

		type = isl_tab_ineq_type(tab, wraps->row[i]);

		if (type == isl_ineq_error)

			return -1;

		if (type == isl_ineq_redundant)

			continue;

		wraps->n_row = 0;

		return 0;

	}

	return 0;

}

/* Return a set that corresponds to the non-redudant constraints

 * (as recorded in tab) of bmap.

 *

 * It's important to remove the redundant constraints as some

 * of the other constraints may have been modified after the

 * constraints were marked redundant.

 * In particular, a constraint may have been relaxed.

 * Redundant constraints are ignored when a constraint is relaxed

 * and should therefore continue to be ignored ever after.

 * Otherwise, the relaxation might be thwarted by some of

 * these constraints.

 */

static __isl_give isl_set *set_from_updated_bmap(__isl_keep isl_basic_map *bmap,

	struct isl_tab *tab)

{

	bmap = isl_basic_map_copy(bmap);

	bmap = isl_basic_map_cow(bmap);

	bmap = isl_basic_map_update_from_tab(bmap, tab);

	return isl_set_from_basic_set(isl_basic_map_underlying_set(bmap));

}

/* Given a basic set i with a constraint k that is adjacent to either the

 * whole of basic set j or a facet of basic set j, check if we can wrap

 * both the facet corresponding to k and the facet of j (or the whole of j)

 * around their ridges to include the other set.

 * If so, replace the pair of basic sets by their union.

 *

 * All constraints of i (except k) are assumed to be valid for j.

 *

 * However, the constraints of j may not be valid for i and so

 * we have to check that the wrapping constraints for j are valid for i.

 *

 * In the case where j has a facet adjacent to i, tab[j] is assumed

 * to have been restricted to this facet, so that the non-redundant

 * constraints in tab[j] are the ridges of the facet.

 * Note that for the purpose of wrapping, it does not matter whether

 * we wrap the ridges of i around the whole of j or just around

 * the facet since all the other constraints are assumed to be valid for j.

 * In practice, we wrap to include the whole of j.

 *        ____			  _____

 *       /    | 		 /     \

 *      /     ||  		/      |

 *      \     ||   	=>	\      |

 *       \    ||		 \     |

 *        \___||		  \____|

 *

 */

static int can_wrap_in_facet(struct isl_map *map, int i, int j, int k,

	struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)

{

	int changed = 0;

	struct isl_wraps wraps;

	isl_mat *mat;

	struct isl_set *set_i = NULL;

	struct isl_set *set_j = NULL;

	struct isl_vec *bound = NULL;

	unsigned total = isl_basic_map_total_dim(map->p[i]);

	struct isl_tab_undo *snap;

	int n;

	set_i = set_from_updated_bmap(map->p[i], tabs[i]);

	set_j = set_from_updated_bmap(map->p[j], tabs[j]);

	mat = isl_mat_alloc(map->ctx, 2 * (map->p[i]->n_eq + map->p[j]->n_eq) +

					map->p[i]->n_ineq + map->p[j]->n_ineq,

					1 + total);

	wraps_init(&wraps, mat, map, i, j, eq_i, ineq_i, eq_j, ineq_j);

	bound = isl_vec_alloc(map->ctx, 1 + total);

	if (!set_i || !set_j || !wraps.mat || !bound)

		goto error;

	isl_seq_cpy(bound->el, map->p[i]->ineq[k], 1 + total);

	isl_int_add_ui(bound->el[0], bound->el[0], 1);

	isl_seq_cpy(wraps.mat->row[0], bound->el, 1 + total);

	wraps.mat->n_row = 1;

	if (add_wraps(&wraps, map->p[j], tabs[j], bound->el, set_i) < 0)

		goto error;

	if (!wraps.mat->n_row)

		goto unbounded;

	snap = isl_tab_snap(tabs[i]);

	if (isl_tab_select_facet(tabs[i], map->p[i]->n_eq + k) < 0)

		goto error;

	if (isl_tab_detect_redundant(tabs[i]) < 0)

		goto error;

	isl_seq_neg(bound->el, map->p[i]->ineq[k], 1 + total);

	n = wraps.mat->n_row;

	if (add_wraps(&wraps, map->p[i], tabs[i], bound->el, set_j) < 0)

		goto error;

	if (isl_tab_rollback(tabs[i], snap) < 0)

		goto error;

	if (check_wraps(wraps.mat, n, tabs[i]) < 0)

		goto error;

	if (!wraps.mat->n_row)

		goto unbounded;

	changed = fuse(map, i, j, tabs, eq_i, ineq_i, eq_j, ineq_j, wraps.mat);

unbounded:

	wraps_free(&wraps);

	isl_set_free(set_i);

	isl_set_free(set_j);

	isl_vec_free(bound);

	return changed;

error:

	wraps_free(&wraps);

	isl_vec_free(bound);

	isl_set_free(set_i);

	isl_set_free(set_j);

	return -1;

}

/* Set the is_redundant property of the "n" constraints in "cuts",

 * except "k" to "v".

 * This is a fairly tricky operation as it bypasses isl_tab.c.

 * The reason we want to temporarily mark some constraints redundant

 * is that we want to ignore them in add_wraps.

 *

 * Initially all cut constraints are non-redundant, but the

 * selection of a facet right before the call to this function

 * may have made some of them redundant.

 * Likewise, the same constraints are marked non-redundant

 * in the second call to this function, before they are officially

 * made non-redundant again in the subsequent rollback.

 */

static void set_is_redundant(struct isl_tab *tab, unsigned n_eq,

	int *cuts, int n, int k, int v)

{

	int l;

	for (l = 0; l < n; ++l) {

		if (l == k)

			continue;

		tab->con[n_eq + cuts[l]].is_redundant = v;

	}

}