/*
 * Copyright 2010      INRIA Saclay
 *
 * Use of this software is governed by the MIT license
 *
 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
 * 91893 Orsay, France 
 */

#include <isl_map_private.h>
#include <isl/set.h>
#include <isl_space_private.h>
#include <isl_seq.h>
#include <isl_aff_private.h>
#include <isl_mat_private.h>
#include <isl_factorization.h>

/*
 * Let C be a cone and define
 *
 *	C' := { y | forall x in C : y x >= 0 }
 *
 * C' contains the coefficients of all linear constraints
 * that are valid for C.
 * Furthermore, C'' = C.
 *
 * If C is defined as { x | A x >= 0 }
 * then any element in C' must be a non-negative combination
 * of the rows of A, i.e., y = t A with t >= 0.  That is,
 *
 *	C' = { y | exists t >= 0 : y = t A }
 *
 * If any of the rows in A actually represents an equality, then
 * also negative combinations of this row are allowed and so the
 * non-negativity constraint on the corresponding element of t
 * can be dropped.
 *
 * A polyhedron P = { x | b + A x >= 0 } can be represented
 * in homogeneous coordinates by the cone
 * C = { [z,x] | b z + A x >= and z >= 0 }
 * The valid linear constraints on C correspond to the valid affine
 * constraints on P.
 * This is essentially Farkas' lemma.
 *
 * Since
 *				  [ 1 0 ]
 *		[ w y ] = [t_0 t] [ b A ]
 *
 * we have
 *
 *	C' = { w, y | exists t_0, t >= 0 : y = t A and w = t_0 + t b }
 * or
 *
 *	C' = { w, y | exists t >= 0 : y = t A and w - t b >= 0 }
 *
 * In practice, we introduce an extra variable (w), shifting all
 * other variables to the right, and an extra inequality
 * (w - t b >= 0) corresponding to the positivity constraint on
 * the homogeneous coordinate.
 *
 * When going back from coefficients to solutions, we immediately
 * plug in 1 for z, which corresponds to shifting all variables
 * to the left, with the leftmost ending up in the constant position.
 */

/* Add the given prefix to all named isl_dim_set dimensions in "space".
 */
static __isl_give isl_space *isl_space_prefix(__isl_take isl_space *space,
	const char *prefix)
{
	int i;
	isl_ctx *ctx;
	isl_size nvar;
	size_t prefix_len = strlen(prefix);

	if (!space)
		return NULL;

	ctx = isl_space_get_ctx(space);
	nvar = isl_space_dim(space, isl_dim_set);
	if (nvar < 0)
		return isl_space_free(space);

	for (i = 0; i < nvar; ++i) {
		const char *name;
		char *prefix_name;

		name = isl_space_get_dim_name(space, isl_dim_set, i);
		if (!name)
			continue;

		prefix_name = isl_alloc_array(ctx, char,
					      prefix_len + strlen(name) + 1);
		if (!prefix_name)
			goto error;
		memcpy(prefix_name, prefix, prefix_len);
		strcpy(prefix_name + prefix_len, name);

		space = isl_space_set_dim_name(space,
						isl_dim_set, i, prefix_name);
		free(prefix_name);
	}

	return space;
error:
	isl_space_free(space);
	return NULL;
}

/* Given a dimension specification of the solutions space, construct
 * a dimension specification for the space of coefficients.
 *
 * In particular transform
 *
 *	[params] -> { S }
 *
 * to
 *
 *	{ coefficients[[cst, params] -> S] }
 *
 * and prefix each dimension name with "c_".
 */
static __isl_give isl_space *isl_space_coefficients(__isl_take isl_space *space)
{
	isl_space *space_param;
	isl_size nvar;
	isl_size nparam;

	nvar = isl_space_dim(space, isl_dim_set);
	nparam = isl_space_dim(space, isl_dim_param);
	if (nvar < 0 || nparam < 0)
		return isl_space_free(space);
	space_param = isl_space_copy(space);
	space_param = isl_space_drop_dims(space_param, isl_dim_set, 0, nvar);
	space_param = isl_space_move_dims(space_param, isl_dim_set, 0,
				 isl_dim_param, 0, nparam);
	space_param = isl_space_prefix(space_param, "c_");
	space_param = isl_space_insert_dims(space_param, isl_dim_set, 0, 1);
	space_param = isl_space_set_dim_name(space_param,
				isl_dim_set, 0, "c_cst");
	space = isl_space_drop_dims(space, isl_dim_param, 0, nparam);
	space = isl_space_prefix(space, "c_");
	space = isl_space_join(isl_space_from_domain(space_param),
			   isl_space_from_range(space));
	space = isl_space_wrap(space);
	space = isl_space_set_tuple_name(space, isl_dim_set, "coefficients");

	return space;
}

/* Drop the given prefix from all named dimensions of type "type" in "space".
 */
static __isl_give isl_space *isl_space_unprefix(__isl_take isl_space *space,
	enum isl_dim_type type, const char *prefix)
{
	int i;
	isl_size n;
	size_t prefix_len = strlen(prefix);

	n = isl_space_dim(space, type);
	if (n < 0)
		return isl_space_free(space);

	for (i = 0; i < n; ++i) {
		const char *name;

		name = isl_space_get_dim_name(space, type, i);
		if (!name)
			continue;
		if (strncmp(name, prefix, prefix_len))
			continue;

		space = isl_space_set_dim_name(space,
						type, i, name + prefix_len);
	}

	return space;
}

/* Given a dimension specification of the space of coefficients, construct
 * a dimension specification for the space of solutions.
 *
 * In particular transform
 *
 *	{ coefficients[[cst, params] -> S] }
 *
 * to
 *
 *	[params] -> { S }
 *
 * and drop the "c_" prefix from the dimension names.
 */
static __isl_give isl_space *isl_space_solutions(__isl_take isl_space *space)
{
	isl_size nparam;

	space = isl_space_unwrap(space);
	space = isl_space_drop_dims(space, isl_dim_in, 0, 1);
	space = isl_space_unprefix(space, isl_dim_in, "c_");
	space = isl_space_unprefix(space, isl_dim_out, "c_");
	nparam = isl_space_dim(space, isl_dim_in);
	if (nparam < 0)
		return isl_space_free(space);
	space = isl_space_move_dims(space,
				    isl_dim_param, 0, isl_dim_in, 0, nparam);
	space = isl_space_range(space);

	return space;
}

/* Return the rational universe basic set in the given space.
 */
static __isl_give isl_basic_set *rational_universe(__isl_take isl_space *space)
{
	isl_basic_set *bset;

	bset = isl_basic_set_universe(space);
	bset = isl_basic_set_set_rational(bset);

	return bset;
}

/* Compute the dual of "bset" by applying Farkas' lemma.
 * As explained above, we add an extra dimension to represent
 * the coefficient of the constant term when going from solutions
 * to coefficients (shift == 1) and we drop the extra dimension when going
 * in the opposite direction (shift == -1).
 * The dual can be created in an arbitrary space.
 * The caller is responsible for putting the result in the appropriate space.
 *
 * If "bset" is (obviously) empty, then the way this emptiness
 * is represented by the constraints does not allow for the application
 * of the standard farkas algorithm.  We therefore handle this case
 * specifically and return the universe basic set.
 */
static __isl_give isl_basic_set *farkas(__isl_take isl_basic_set *bset,
	int shift)
{
	int i, j, k;
	isl_ctx *ctx;
	isl_space *space;
	isl_basic_set *dual = NULL;
	isl_size total;

	total = isl_basic_set_dim(bset, isl_dim_all);
	if (total < 0)
		return isl_basic_set_free(bset);

	ctx = isl_basic_set_get_ctx(bset);
	space = isl_space_set_alloc(ctx, 0, total + shift);
	if (isl_basic_set_plain_is_empty(bset)) {
		isl_basic_set_free(bset);
		return rational_universe(space);
	}

	dual = isl_basic_set_alloc_space(space, bset->n_eq + bset->n_ineq,
					total, bset->n_ineq + (shift > 0));
	dual = isl_basic_set_set_rational(dual);

	for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
		k = isl_basic_set_alloc_div(dual);
		if (k < 0)
			goto error;
		isl_int_set_si(dual->div[k][0], 0);
	}

	for (i = 0; i < total; ++i) {
		k = isl_basic_set_alloc_equality(dual);
		if (k < 0)
			goto error;
		isl_seq_clr(dual->eq[k], 1 + shift + total);
		isl_int_set_si(dual->eq[k][1 + shift + i], -1);
		for (j = 0; j < bset->n_eq; ++j)
			isl_int_set(dual->eq[k][1 + shift + total + j],
				    bset->eq[j][1 + i]);
		for (j = 0; j < bset->n_ineq; ++j)
			isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
				    bset->ineq[j][1 + i]);
	}

	for (i = 0; i < bset->n_ineq; ++i) {
		k = isl_basic_set_alloc_inequality(dual);
		if (k < 0)
			goto error;
		isl_seq_clr(dual->ineq[k],
			    1 + shift + total + bset->n_eq + bset->n_ineq);
		isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
	}

	if (shift > 0) {
		k = isl_basic_set_alloc_inequality(dual);
		if (k < 0)
			goto error;
		isl_seq_clr(dual->ineq[k], 2 + total);
		isl_int_set_si(dual->ineq[k][1], 1);
		for (j = 0; j < bset->n_eq; ++j)
			isl_int_neg(dual->ineq[k][2 + total + j],
				    bset->eq[j][0]);
		for (j = 0; j < bset->n_ineq; ++j)
			isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
				    bset->ineq[j][0]);
	}

	dual = isl_basic_set_remove_divs(dual);
	dual = isl_basic_set_simplify(dual);
	dual = isl_basic_set_finalize(dual);

	isl_basic_set_free(bset);
	return dual;
error:
	isl_basic_set_free(bset);
	isl_basic_set_free(dual);
	return NULL;
}

/* Construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the given basic set, ignoring
 * the space of input and output and without any further decomposition.
 */
static __isl_give isl_basic_set *isl_basic_set_coefficients_base(
	__isl_take isl_basic_set *bset)
{
	return farkas(bset, 1);
}

/* Return the inverse mapping of "morph".
 */
static __isl_give isl_mat *peek_inv(__isl_keep isl_morph *morph)
{
	return morph ? morph->inv : NULL;
}

/* Return a copy of the inverse mapping of "morph".
 */
static __isl_give isl_mat *get_inv(__isl_keep isl_morph *morph)
{
	return isl_mat_copy(peek_inv(morph));
}

/* Information about a single factor within isl_basic_set_coefficients_product.
 *
 * "start" is the position of the first coefficient (beyond
 * the one corresponding to the constant term) in this factor.
 * "dim" is the number of coefficients (other than
 * the one corresponding to the constant term) in this factor.
 * "n_line" is the number of lines in "coeff".
 * "n_ray" is the number of rays (other than lines) in "coeff".
 * "n_vertex" is the number of vertices in "coeff".
 *
 * While iterating over the vertices,
 * "pos" represents the inequality constraint corresponding
 * to the current vertex.
 */
struct isl_coefficients_factor_data {
	isl_basic_set *coeff;
	int start;
	int dim;
	int n_line;
	int n_ray;
	int n_vertex;
	int pos;
};

/* Internal data structure for isl_basic_set_coefficients_product.
 * "n" is the number of factors in the factorization.
 * "pos" is the next factor that will be considered.
 * "start_next" is the position of the first coefficient (beyond
 * the one corresponding to the constant term) in the next factor.
 * "factors" contains information about the individual "n" factors.
 */
struct isl_coefficients_product_data {
	int n;
	int pos;
	int start_next;
	struct isl_coefficients_factor_data *factors;
};

/* Initialize the internal data structure for
 * isl_basic_set_coefficients_product.
 */
static isl_stat isl_coefficients_product_data_init(isl_ctx *ctx,
	struct isl_coefficients_product_data *data, int n)
{
	data->n = n;
	data->pos = 0;
	data->start_next = 0;
	data->factors = isl_calloc_array(ctx,
					struct isl_coefficients_factor_data, n);
	if (!data->factors)
		return isl_stat_error;
	return isl_stat_ok;
}

/* Free all memory allocated in "data".
 */
static void isl_coefficients_product_data_clear(
	struct isl_coefficients_product_data *data)
{
	int i;

	if (data->factors) {
		for (i = 0; i < data->n; ++i) {
			isl_basic_set_free(data->factors[i].coeff);
		}
	}
	free(data->factors);
}

/* Does inequality "ineq" in the (dual) basic set "bset" represent a ray?
 * In particular, does it have a zero denominator
 * (i.e., a zero coefficient for the constant term)?
 */
static int is_ray(__isl_keep isl_basic_set *bset, int ineq)
{
	return isl_int_is_zero(bset->ineq[ineq][1]);
}

/* isl_factorizer_every_factor_basic_set callback that
 * constructs a basic set containing the tuples of coefficients of all
 * valid affine constraints on the factor "bset" and
 * extracts further information that will be used
 * when combining the results over the different factors.
 */
static isl_bool isl_basic_set_coefficients_factor(
	__isl_keep isl_basic_set *bset, void *user)
{
	struct isl_coefficients_product_data *data = user;
	isl_basic_set *coeff;
	isl_size n_eq, n_ineq, dim;
	int i, n_ray, n_vertex;

	coeff = isl_basic_set_coefficients_base(isl_basic_set_copy(bset));
	data->factors[data->pos].coeff = coeff;
	if (!coeff)
		return isl_bool_error;

	dim = isl_basic_set_dim(bset, isl_dim_set);
	n_eq = isl_basic_set_n_equality(coeff);
	n_ineq = isl_basic_set_n_inequality(coeff);
	if (dim < 0 || n_eq < 0 || n_ineq < 0)
		return isl_bool_error;
	n_ray = n_vertex = 0;
	for (i = 0; i < n_ineq; ++i) {
		if (is_ray(coeff, i))
			n_ray++;
		else
			n_vertex++;
	}
	data->factors[data->pos].start = data->start_next;
	data->factors[data->pos].dim = dim;
	data->factors[data->pos].n_line = n_eq;
	data->factors[data->pos].n_ray = n_ray;
	data->factors[data->pos].n_vertex = n_vertex;
	data->pos++;
	data->start_next += dim;

	return isl_bool_true;
}

/* Clear an entry in the product, given that there is a "total" number
 * of coefficients (other than that of the constant term).
 */
static void clear_entry(isl_int *entry, int total)
{
	isl_seq_clr(entry, 1 + 1 + total);
}

/* Set the part of the entry corresponding to factor "data",
 * from the factor coefficients in "src".
 */
static void set_factor(isl_int *entry, isl_int *src,
	struct isl_coefficients_factor_data *data)
{
	isl_seq_cpy(entry + 1 + 1 + data->start, src + 1 + 1, data->dim);
}

/* Set the part of the entry corresponding to factor "data",
 * from the factor coefficients in "src" multiplied by "f".
 */
static void scale_factor(isl_int *entry, isl_int *src, isl_int f,
	struct isl_coefficients_factor_data *data)
{
	isl_seq_scale(entry + 1 + 1 + data->start, src + 1 + 1, f, data->dim);
}

/* Add all lines from the given factor to "bset",
 * given that there is a "total" number of coefficients
 * (other than that of the constant term).
 */
static __isl_give isl_basic_set *add_lines(__isl_take isl_basic_set *bset,
	struct isl_coefficients_factor_data *factor, int total)
{
	int i;

	for (i = 0; i < factor->n_line; ++i) {
		int k;

		k = isl_basic_set_alloc_equality(bset);
		if (k < 0)
			return isl_basic_set_free(bset);
		clear_entry(bset->eq[k], total);
		set_factor(bset->eq[k], factor->coeff->eq[i], factor);
	}

	return bset;
}

/* Add all rays (other than lines) from the given factor to "bset",
 * given that there is a "total" number of coefficients
 * (other than that of the constant term).
 */
static __isl_give isl_basic_set *add_rays(__isl_take isl_basic_set *bset,
	struct isl_coefficients_factor_data *data, int total)
{
	int i;
	int n_ineq = data->n_ray + data->n_vertex;

	for (i = 0; i < n_ineq; ++i) {
		int k;

		if (!is_ray(data->coeff, i))
			continue;

		k = isl_basic_set_alloc_inequality(bset);
		if (k < 0)
			return isl_basic_set_free(bset);
		clear_entry(bset->ineq[k], total);
		set_factor(bset->ineq[k], data->coeff->ineq[i], data);
	}

	return bset;
}

/* Move to the first vertex of the given factor starting
 * at inequality constraint "start", setting factor->pos and
 * returning 1 if a vertex is found.
 */
static int factor_first_vertex(struct isl_coefficients_factor_data *factor,
	int start)
{
	int j;
	int n = factor->n_ray + factor->n_vertex;

	for (j = start; j < n; ++j) {
		if (is_ray(factor->coeff, j))
			continue;
		factor->pos = j;
		return 1;
	}

	return 0;
}

/* Move to the first constraint in each factor starting at "first"
 * that represents a vertex.
 * In particular, skip the initial constraints that correspond to rays.
 */
static void first_vertex(struct isl_coefficients_product_data *data, int first)
{
	int i;

	for (i = first; i < data->n; ++i)
		factor_first_vertex(&data->factors[i], 0);
}

/* Move to the next vertex in the product.
 * In particular, move to the next vertex of the last factor.
 * If all vertices of this last factor have already been considered,
 * then move to the next vertex of the previous factor(s)
 * until a factor is found that still has a next vertex.
 * Once such a next vertex has been found, the subsequent
 * factors are reset to the first vertex.
 * Return 1 if any next vertex was found.
 */
static int next_vertex(struct isl_coefficients_product_data *data)
{
	int i;

	for (i = data->n - 1; i >= 0; --i) {
		struct isl_coefficients_factor_data *factor = &data->factors[i];

		if (!factor_first_vertex(factor, factor->pos + 1))
			continue;
		first_vertex(data, i + 1);
		return 1;
	}

	return 0;
}

/* Add a vertex to the product "bset" combining the currently selected
 * vertices of the factors.
 *
 * In the dual representation, the constant term is always zero.
 * The vertex itself is the sum of the contributions of the factors
 * with a shared denominator in position 1.
 *
 * First compute the shared denominator (lcm) and
 * then scale the numerators to this shared denominator.
 */
static __isl_give isl_basic_set *add_vertex(__isl_take isl_basic_set *bset,
	struct isl_coefficients_product_data *data)
{
	int i;
	int k;
	isl_int lcm, f;

	k = isl_basic_set_alloc_inequality(bset);
	if (k < 0)
		return isl_basic_set_free(bset);

	isl_int_init(lcm);
	isl_int_init(f);
	isl_int_set_si(lcm, 1);
	for (i = 0; i < data->n; ++i) {
		struct isl_coefficients_factor_data *factor = &data->factors[i];
		isl_basic_set *coeff = factor->coeff;
		int pos = factor->pos;
		isl_int_lcm(lcm, lcm, coeff->ineq[pos][1]);
	}
	isl_int_set_si(bset->ineq[k][0], 0);
	isl_int_set(bset->ineq[k][1], lcm);

	for (i = 0; i < data->n; ++i) {
		struct isl_coefficients_factor_data *factor = &data->factors[i];
		isl_basic_set *coeff = factor->coeff;
		int pos = factor->pos;
		isl_int_divexact(f, lcm, coeff->ineq[pos][1]);
		scale_factor(bset->ineq[k], coeff->ineq[pos], f, factor);
	}

	isl_int_clear(f);
	isl_int_clear(lcm);

	return bset;
}

/* Combine the duals of the factors in the factorization of a basic set
 * to form the dual of the entire basic set.
 * The dual share the coefficient of the constant term.
 * All other coefficients are specific to a factor.
 * Any constraint not involving the coefficient of the constant term
 * can therefor simply be copied into the appropriate position.
 * This includes all equality constraints since the coefficient
 * of the constant term can always be increased and therefore
 * never appears in an equality constraint.
 * The inequality constraints involving the coefficient of
 * the constant term need to be combined across factors.
 * In particular, if this coefficient needs to be greater than or equal
 * to some linear combination of the other coefficients in each factor,
 * then it needs to be greater than or equal to the sum of
 * these linear combinations across the factors.
 *
 * Alternatively, the constraints of the dual can be seen
 * as the vertices, rays and lines of the original basic set.
 * Clearly, rays and lines can simply be copied,
 * while vertices needs to be combined across factors.
 * This means that the number of rays and lines in the product
 * is equal to the sum of the numbers in the factors,
 * while the number of vertices is the product
 * of the number of vertices in the factors.  Note that each
 * factor has at least one vertex.
 * The only exception is when the factor is the dual of an obviously empty set,
 * in which case a universe dual is created.
 * In this case, return a universe dual for the product as well.
 *
 * While constructing the vertices, look for the first combination
 * of inequality constraints that represent a vertex,
 * construct the corresponding vertex and then move on
 * to the next combination of inequality constraints until
 * all combinations have been considered.
 */
static __isl_give isl_basic_set *construct_product(isl_ctx *ctx,
	struct isl_coefficients_product_data *data)
{
	int i;
	int n_line, n_ray, n_vertex;
	int total;
	isl_space *space;
	isl_basic_set *product;

	if (!data->factors)
		return NULL;

	total = data->start_next;

	n_line = 0;
	n_ray = 0;
	n_vertex = 1;
	for (i = 0; i < data->n; ++i) {
		n_line += data->factors[i].n_line;
		n_ray += data->factors[i].n_ray;
		n_vertex *= data->factors[i].n_vertex;
	}

	space = isl_space_set_alloc(ctx, 0, 1 + total);
	if (n_vertex == 0)
		return rational_universe(space);
	product = isl_basic_set_alloc_space(space, 0, n_line, n_ray + n_vertex);
	product = isl_basic_set_set_rational(product);

	for (i = 0; i < data->n; ++i)
		product = add_lines(product, &data->factors[i], total);
	for (i = 0; i < data->n; ++i)
		product = add_rays(product, &data->factors[i], total);

	first_vertex(data, 0);
	do {
		product = add_vertex(product, data);
	} while (next_vertex(data));

	return product;
}

/* Given a factorization "f" of a basic set,
 * construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the product of the factors, ignoring
 * the space of input and output.
 * Note that this product may not be equal to the original basic set,
 * if a non-trivial transformation is involved.
 * This is handled by the caller.
 *
 * Compute the tuples of coefficients for each factor separately and
 * then combine the results.
 */
static __isl_give isl_basic_set *isl_basic_set_coefficients_product(
	__isl_take isl_factorizer *f)
{
	struct isl_coefficients_product_data data;
	isl_ctx *ctx;
	isl_basic_set *coeff;
	isl_bool every;

	ctx = isl_factorizer_get_ctx(f);
	if (isl_coefficients_product_data_init(ctx, &data, f->n_group) < 0)
		f = isl_factorizer_free(f);
	every = isl_factorizer_every_factor_basic_set(f,
			&isl_basic_set_coefficients_factor, &data);
	isl_factorizer_free(f);
	if (every >= 0)
		coeff = construct_product(ctx, &data);
	else
		coeff = NULL;
	isl_coefficients_product_data_clear(&data);

	return coeff;
}

/* Given a factorization "f" of a basic set,
 * construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the basic set, ignoring
 * the space of input and output.
 *
 * The factorization may involve a linear transformation of the basic set.
 * In particular, the transformed basic set is formulated
 * in terms of x' = U x, i.e., x = V x', with V = U^{-1}.
 * The dual is then computed in terms of y' with y'^t [z; x'] >= 0.
 * Plugging in y' = [1 0; 0 V^t] y yields
 * y^t [1 0; 0 V] [z; x'] >= 0, i.e., y^t [z; x] >= 0, which is
 * the desired set of coefficients y.
 * Note that this transformation to y' only needs to be applied
 * if U is not the identity matrix.
 */
static __isl_give isl_basic_set *isl_basic_set_coefficients_morphed_product(
	__isl_take isl_factorizer *f)
{
	isl_bool is_identity;
	isl_space *space;
	isl_mat *inv;
	isl_multi_aff *ma;
	isl_basic_set *coeff;

	if (!f)
		goto error;
	is_identity = isl_mat_is_scaled_identity(peek_inv(f->morph));
	if (is_identity < 0)
		goto error;
	if (is_identity)
		return isl_basic_set_coefficients_product(f);

	inv = get_inv(f->morph);
	inv = isl_mat_transpose(inv);
	inv = isl_mat_lin_to_aff(inv);

	coeff = isl_basic_set_coefficients_product(f);
	space = isl_space_map_from_set(isl_basic_set_get_space(coeff));
	ma = isl_multi_aff_from_aff_mat(space, inv);
	coeff = isl_basic_set_preimage_multi_aff(coeff, ma);

	return coeff;
error:
	isl_factorizer_free(f);
	return NULL;
}

/* Construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the given basic set, ignoring
 * the space of input and output.
 *
 * The caller has already checked that "bset" does not involve
 * any local variables.  It may have parameters, though.
 * Treat them as regular variables internally.
 * This is especially important for the factorization,
 * since the (original) parameters should be taken into account
 * explicitly in this factorization.
 *
 * Check if the basic set can be factorized.
 * If so, compute constraints on the coefficients of the factors
 * separately and combine the results.
 * Otherwise, compute the results for the input basic set as a whole.
 */
static __isl_give isl_basic_set *basic_set_coefficients(
	__isl_take isl_basic_set *bset)
{
	isl_factorizer *f;
	isl_size nparam;

	nparam = isl_basic_set_dim(bset, isl_dim_param);
	if (nparam < 0)
		return isl_basic_set_free(bset);
	bset = isl_basic_set_move_dims(bset, isl_dim_set, 0,
					    isl_dim_param, 0, nparam);

	f = isl_basic_set_factorizer(bset);
	if (!f)
		return isl_basic_set_free(bset);
	if (f->n_group > 0) {
		isl_basic_set_free(bset);
		return isl_basic_set_coefficients_morphed_product(f);
	}
	isl_factorizer_free(f);
	return isl_basic_set_coefficients_base(bset);
}

/* Construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the given basic set.
 */
__isl_give isl_basic_set *isl_basic_set_coefficients(
	__isl_take isl_basic_set *bset)
{
	isl_space *space;

	if (!bset)
		return NULL;
	if (bset->n_div)
		isl_die(bset->ctx, isl_error_invalid,
			"input set not allowed to have local variables",
			goto error);

	space = isl_basic_set_get_space(bset);
	space = isl_space_coefficients(space);

	bset = basic_set_coefficients(bset);
	bset = isl_basic_set_reset_space(bset, space);
	return bset;
error:
	isl_basic_set_free(bset);
	return NULL;
}

/* Construct a basic set containing the elements that satisfy all
 * affine constraints whose coefficient tuples are
 * contained in the given basic set.
 */
__isl_give isl_basic_set *isl_basic_set_solutions(
	__isl_take isl_basic_set *bset)
{
	isl_space *space;

	if (!bset)
		return NULL;
	if (bset->n_div)
		isl_die(bset->ctx, isl_error_invalid,
			"input set not allowed to have local variables",
			goto error);

	space = isl_basic_set_get_space(bset);
	space = isl_space_solutions(space);

	bset = farkas(bset, -1);
	bset = isl_basic_set_reset_space(bset, space);
	return bset;
error:
	isl_basic_set_free(bset);
	return NULL;
}

/* Construct a basic set containing the tuples of coefficients of all
 * valid affine constraints on the given set.
 */
__isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
{
	int i;
	isl_basic_set *coeff;

	if (!set)
		return NULL;
	if (set->n == 0) {
		isl_space *space = isl_set_get_space(set);
		space = isl_space_coefficients(space);
		isl_set_free(set);
		return rational_universe(space);
	}

	coeff = isl_basic_set_coefficients(isl_basic_set_copy(set->p[0]));

	for (i = 1; i < set->n; ++i) {
		isl_basic_set *bset, *coeff_i;
		bset = isl_basic_set_copy(set->p[i]);
		coeff_i = isl_basic_set_coefficients(bset);
		coeff = isl_basic_set_intersect(coeff, coeff_i);
	}

	isl_set_free(set);
	return coeff;
}

/* Wrapper around isl_basic_set_coefficients for use
 * as a isl_basic_set_list_map callback.
 */
static __isl_give isl_basic_set *coefficients_wrap(
	__isl_take isl_basic_set *bset, void *user)
{
	return isl_basic_set_coefficients(bset);
}

/* Replace the elements of "list" by the result of applying
 * isl_basic_set_coefficients to them.
 */
__isl_give isl_basic_set_list *isl_basic_set_list_coefficients(
	__isl_take isl_basic_set_list *list)
{
	return isl_basic_set_list_map(list, &coefficients_wrap, NULL);
}

/* Construct a basic set containing the elements that satisfy all
 * affine constraints whose coefficient tuples are
 * contained in the given set.
 */
__isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
{
	int i;
	isl_basic_set *sol;

	if (!set)
		return NULL;
	if (set->n == 0) {
		isl_space *space = isl_set_get_space(set);
		space = isl_space_solutions(space);
		isl_set_free(set);
		return rational_universe(space);
	}

	sol = isl_basic_set_solutions(isl_basic_set_copy(set->p[0]));

	for (i = 1; i < set->n; ++i) {
		isl_basic_set *bset, *sol_i;
		bset = isl_basic_set_copy(set->p[i]);
		sol_i = isl_basic_set_solutions(bset);
		sol = isl_basic_set_intersect(sol, sol_i);
	}

	isl_set_free(set);
	return sol;
}