/*

 * 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 

 */

#define ISL_DIM_H

#include <isl_map_private.h>

#include <isl_union_map_private.h>

#include <isl_polynomial_private.h>

#include <isl_point_private.h>

#include <isl_space_private.h>

#include <isl_lp_private.h>

#include <isl_seq.h>

#include <isl_mat_private.h>

#include <isl_val_private.h>

#include <isl_vec_private.h>

#include <isl_config.h>

#include <isl/deprecated/polynomial_int.h>

enum isl_fold isl_fold_type_negate(enum isl_fold type)

{

	switch (type) {

	case isl_fold_min:

		return isl_fold_max;

	case isl_fold_max:

		return isl_fold_min;

	case isl_fold_list:

		return isl_fold_list;

	}

	isl_die(NULL, isl_error_internal, "unhandled isl_fold type", abort());

}

static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(

	enum isl_fold type, __isl_take isl_space *dim, int n)

{

	isl_qpolynomial_fold *fold;

	if (!dim)

		goto error;

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

	fold = isl_calloc(dim->ctx, struct isl_qpolynomial_fold,

			sizeof(struct isl_qpolynomial_fold) +

			(n - 1) * sizeof(struct isl_qpolynomial *));

	if (!fold)

		goto error;

	fold->ref = 1;

	fold->size = n;

	fold->n = 0;

	fold->type = type;

	fold->dim = dim;

	return fold;

error:

	isl_space_free(dim);

	return NULL;

}

isl_ctx *isl_qpolynomial_fold_get_ctx(__isl_keep isl_qpolynomial_fold *fold)

{

	return fold ? fold->dim->ctx : NULL;

}

__isl_give isl_space *isl_qpolynomial_fold_get_domain_space(

	__isl_keep isl_qpolynomial_fold *fold)

{

	return fold ? isl_space_copy(fold->dim) : NULL;

}

__isl_give isl_space *isl_qpolynomial_fold_get_space(

	__isl_keep isl_qpolynomial_fold *fold)

{

	isl_space *space;

	if (!fold)

		return NULL;

	space = isl_space_copy(fold->dim);

	space = isl_space_from_domain(space);

	space = isl_space_add_dims(space, isl_dim_out, 1);

	return space;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_reset_domain_space(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_space *dim)

{

	int i;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold || !dim)

		goto error;

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

		fold->qp[i] = isl_qpolynomial_reset_domain_space(fold->qp[i],

							isl_space_copy(dim));

		if (!fold->qp[i])

			goto error;

	}

	isl_space_free(fold->dim);

	fold->dim = dim;

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	isl_space_free(dim);

	return NULL;

}

/* Reset the space of "fold".  This function is called from isl_pw_templ.c

 * and doesn't know if the space of an element object is represented

 * directly or through its domain.  It therefore passes along both.

 */

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_reset_space_and_domain(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_space *space,

	__isl_take isl_space *domain)

{

	isl_space_free(space);

	return isl_qpolynomial_fold_reset_domain_space(fold, domain);

}

int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,

	enum isl_dim_type type, unsigned first, unsigned n)

{

	int i;

	if (!fold)

		return -1;

	if (fold->n == 0 || n == 0)

		return 0;

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

		int involves = isl_qpolynomial_involves_dims(fold->qp[i],

							    type, first, n);

		if (involves < 0 || involves)

			return involves;

	}

	return 0;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_set_dim_name(

	__isl_take isl_qpolynomial_fold *fold,

	enum isl_dim_type type, unsigned pos, const char *s)

{

	int i;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		return NULL;

	fold->dim = isl_space_set_dim_name(fold->dim, type, pos, s);

	if (!fold->dim)

		goto error;

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

		fold->qp[i] = isl_qpolynomial_set_dim_name(fold->qp[i],

							    type, pos, s);

		if (!fold->qp[i])

			goto error;

	}

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(

	__isl_take isl_qpolynomial_fold *fold,

	enum isl_dim_type type, unsigned first, unsigned n)

{

	int i;

	enum isl_dim_type set_type;

	if (!fold)

		return NULL;

	if (n == 0)

		return fold;

	set_type = type == isl_dim_in ? isl_dim_set : type;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		return NULL;

	fold->dim = isl_space_drop_dims(fold->dim, set_type, first, n);

	if (!fold->dim)

		goto error;

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

		fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],

							    type, first, n);

		if (!fold->qp[i])

			goto error;

	}

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_insert_dims(

	__isl_take isl_qpolynomial_fold *fold,

	enum isl_dim_type type, unsigned first, unsigned n)

{

	int i;

	if (!fold)

		return NULL;

	if (n == 0 && !isl_space_is_named_or_nested(fold->dim, type))

		return fold;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		return NULL;

	fold->dim = isl_space_insert_dims(fold->dim, type, first, n);

	if (!fold->dim)

		goto error;

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

		fold->qp[i] = isl_qpolynomial_insert_dims(fold->qp[i],

							    type, first, n);

		if (!fold->qp[i])

			goto error;

	}

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	return NULL;

}

static int isl_qpolynomial_cst_sign(__isl_keep isl_qpolynomial *qp)

{

	struct isl_upoly_cst *cst;

	cst = isl_upoly_as_cst(qp->upoly);

	if (!cst)

		return 0;

	return isl_int_sgn(cst->n) < 0 ? -1 : 1;

}

static int isl_qpolynomial_aff_sign(__isl_keep isl_set *set,

	__isl_keep isl_qpolynomial *qp)

{

	enum isl_lp_result res;

	isl_vec *aff;

	isl_int opt;

	int sgn = 0;

	aff = isl_qpolynomial_extract_affine(qp);

	if (!aff)

		return 0;

	isl_int_init(opt);

	res = isl_set_solve_lp(set, 0, aff->el + 1, aff->el[0],

				&opt, NULL, NULL);

	if (res == isl_lp_error)

		goto done;

	if (res == isl_lp_empty ||

	    (res == isl_lp_ok && !isl_int_is_neg(opt))) {

		sgn = 1;

		goto done;

	}

	res = isl_set_solve_lp(set, 1, aff->el + 1, aff->el[0],

				&opt, NULL, NULL);

	if (res == isl_lp_ok && !isl_int_is_pos(opt))

		sgn = -1;

done:

	isl_int_clear(opt);

	isl_vec_free(aff);

	return sgn;

}

/* Determine, if possible, the sign of the quasipolynomial "qp" on

 * the domain "set".

 *

 * If qp is a constant, then the problem is trivial.

 * If qp is linear, then we check if the minimum of the corresponding

 * affine constraint is non-negative or if the maximum is non-positive.

 *

 * Otherwise, we check if the outermost variable "v" has a lower bound "l"

 * in "set".  If so, we write qp(v,v') as

 *

 *	q(v,v') * (v - l) + r(v')

 *

 * if q(v,v') and r(v') have the same known sign, then the original

 * quasipolynomial has the same sign as well.

 *

 * Return

 *	-1 if qp <= 0

 *	 1 if qp >= 0

 *	 0 if unknown

 */

static int isl_qpolynomial_sign(__isl_keep isl_set *set,

	__isl_keep isl_qpolynomial *qp)

{

	int d;

	int i;

	int is;

	struct isl_upoly_rec *rec;

	isl_vec *v;

	isl_int l;

	enum isl_lp_result res;

	int sgn = 0;

	is = isl_qpolynomial_is_cst(qp, NULL, NULL);

	if (is < 0)

		return 0;

	if (is)

		return isl_qpolynomial_cst_sign(qp);

	is = isl_qpolynomial_is_affine(qp);

	if (is < 0)

		return 0;

	if (is)

		return isl_qpolynomial_aff_sign(set, qp);

	if (qp->div->n_row > 0)

		return 0;

	rec = isl_upoly_as_rec(qp->upoly);

	if (!rec)

		return 0;

	d = isl_space_dim(qp->dim, isl_dim_all);

	v = isl_vec_alloc(set->ctx, 2 + d);

	if (!v)

		return 0;

	isl_seq_clr(v->el + 1, 1 + d);

	isl_int_set_si(v->el[0], 1);

	isl_int_set_si(v->el[2 + qp->upoly->var], 1);

	isl_int_init(l);

	res = isl_set_solve_lp(set, 0, v->el + 1, v->el[0], &l, NULL, NULL);

	if (res == isl_lp_ok) {

		isl_qpolynomial *min;

		isl_qpolynomial *base;

		isl_qpolynomial *r, *q;

		isl_qpolynomial *t;

		min = isl_qpolynomial_cst_on_domain(isl_space_copy(qp->dim), l);

		base = isl_qpolynomial_var_pow_on_domain(isl_space_copy(qp->dim),

						qp->upoly->var, 1);

		r = isl_qpolynomial_alloc(isl_space_copy(qp->dim), 0,

					  isl_upoly_copy(rec->p[rec->n - 1]));

		q = isl_qpolynomial_copy(r);

		for (i = rec->n - 2; i >= 0; --i) {

			r = isl_qpolynomial_mul(r, isl_qpolynomial_copy(min));

			t = isl_qpolynomial_alloc(isl_space_copy(qp->dim), 0,

						  isl_upoly_copy(rec->p[i]));

			r = isl_qpolynomial_add(r, t);

			if (i == 0)

				break;

			q = isl_qpolynomial_mul(q, isl_qpolynomial_copy(base));

			q = isl_qpolynomial_add(q, isl_qpolynomial_copy(r));

		}

		if (isl_qpolynomial_is_zero(q))

			sgn = isl_qpolynomial_sign(set, r);

		else if (isl_qpolynomial_is_zero(r))

			sgn = isl_qpolynomial_sign(set, q);

		else {

			int sgn_q, sgn_r;

			sgn_r = isl_qpolynomial_sign(set, r);

			sgn_q = isl_qpolynomial_sign(set, q);

			if (sgn_r == sgn_q)

				sgn = sgn_r;

		}

		isl_qpolynomial_free(min);

		isl_qpolynomial_free(base);

		isl_qpolynomial_free(q);

		isl_qpolynomial_free(r);

	}

	isl_int_clear(l);

	isl_vec_free(v);

	return sgn;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold_on_domain(

	__isl_keep isl_set *set,

	__isl_take isl_qpolynomial_fold *fold1,

	__isl_take isl_qpolynomial_fold *fold2)

{

	int i, j;

	int n1;

	struct isl_qpolynomial_fold *res = NULL;

	int better;

	if (!fold1 || !fold2)

		goto error;

	isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);

	isl_assert(fold1->dim->ctx, isl_space_is_equal(fold1->dim, fold2->dim),

			goto error);

	better = fold1->type == isl_fold_max ? -1 : 1;

	if (isl_qpolynomial_fold_is_empty(fold1)) {

		isl_qpolynomial_fold_free(fold1);

		return fold2;

	}

	if (isl_qpolynomial_fold_is_empty(fold2)) {

		isl_qpolynomial_fold_free(fold2);

		return fold1;

	}

	res = qpolynomial_fold_alloc(fold1->type, isl_space_copy(fold1->dim),

					fold1->n + fold2->n);

	if (!res)

		goto error;

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

		res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);

		if (!res->qp[res->n])

			goto error;

		res->n++;

	}

	n1 = res->n;

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

		for (j = n1 - 1; j >= 0; --j) {

			isl_qpolynomial *d;

			int sgn;

			d = isl_qpolynomial_sub(

				isl_qpolynomial_copy(res->qp[j]),

				isl_qpolynomial_copy(fold2->qp[i]));

			sgn = isl_qpolynomial_sign(set, d);

			isl_qpolynomial_free(d);

			if (sgn == 0)

				continue;

			if (sgn != better)

				break;

			isl_qpolynomial_free(res->qp[j]);

			if (j != n1 - 1)

				res->qp[j] = res->qp[n1 - 1];

			n1--;

			if (n1 != res->n - 1)

				res->qp[n1] = res->qp[res->n - 1];

			res->n--;

		}

		if (j >= 0)

			continue;

		res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);

		if (!res->qp[res->n])

			goto error;

		res->n++;

	}

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return res;

error:

	isl_qpolynomial_fold_free(res);

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_add_qpolynomial(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_qpolynomial *qp)

{

	int i;

	if (!fold || !qp)

		goto error;

	if (isl_qpolynomial_is_zero(qp)) {

		isl_qpolynomial_free(qp);

		return fold;

	}

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		goto error;

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

		fold->qp[i] = isl_qpolynomial_add(fold->qp[i],

						isl_qpolynomial_copy(qp));

		if (!fold->qp[i])

			goto error;

	}

	isl_qpolynomial_free(qp);

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	isl_qpolynomial_free(qp);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_add_on_domain(

	__isl_keep isl_set *dom,

	__isl_take isl_qpolynomial_fold *fold1,

	__isl_take isl_qpolynomial_fold *fold2)

{

	int i;

	isl_qpolynomial_fold *res = NULL;

	if (!fold1 || !fold2)

		goto error;

	if (isl_qpolynomial_fold_is_empty(fold1)) {

		isl_qpolynomial_fold_free(fold1);

		return fold2;

	}

	if (isl_qpolynomial_fold_is_empty(fold2)) {

		isl_qpolynomial_fold_free(fold2);

		return fold1;

	}

	if (fold1->n == 1 && fold2->n != 1)

		return isl_qpolynomial_fold_add_on_domain(dom, fold2, fold1);

	if (fold2->n == 1) {

		res = isl_qpolynomial_fold_add_qpolynomial(fold1,

					isl_qpolynomial_copy(fold2->qp[0]));

		isl_qpolynomial_fold_free(fold2);

		return res;

	}

	res = isl_qpolynomial_fold_add_qpolynomial(

				isl_qpolynomial_fold_copy(fold1),

				isl_qpolynomial_copy(fold2->qp[0]));

	for (i = 1; i < fold2->n; ++i) {

		isl_qpolynomial_fold *res_i;

		res_i = isl_qpolynomial_fold_add_qpolynomial(

					isl_qpolynomial_fold_copy(fold1),

					isl_qpolynomial_copy(fold2->qp[i]));

		res = isl_qpolynomial_fold_fold_on_domain(dom, res, res_i);

	}

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return res;

error:

	isl_qpolynomial_fold_free(res);

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_substitute_equalities(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_basic_set *eq)

{

	int i;

	if (!fold || !eq)

		goto error;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		return NULL;

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

		fold->qp[i] = isl_qpolynomial_substitute_equalities(fold->qp[i],

							isl_basic_set_copy(eq));

		if (!fold->qp[i])

			goto error;

	}

	isl_basic_set_free(eq);

	return fold;

error:

	isl_basic_set_free(eq);

	isl_qpolynomial_fold_free(fold);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_gist(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *context)

{

	int i;

	if (!fold || !context)

		goto error;

	fold = isl_qpolynomial_fold_cow(fold);

	if (!fold)

		return NULL;

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

		fold->qp[i] = isl_qpolynomial_gist(fold->qp[i],

							isl_set_copy(context));

		if (!fold->qp[i])

			goto error;

	}

	isl_set_free(context);

	return fold;

error:

	isl_set_free(context);

	isl_qpolynomial_fold_free(fold);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_gist_params(

	__isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *context)

{

	isl_space *space = isl_qpolynomial_fold_get_domain_space(fold);

	isl_set *dom_context = isl_set_universe(space);

	dom_context = isl_set_intersect_params(dom_context, context);

	return isl_qpolynomial_fold_gist(fold, dom_context);

}

#define HAS_TYPE

#undef PW

#define PW isl_pw_qpolynomial_fold

#undef EL

#define EL isl_qpolynomial_fold

#undef EL_IS_ZERO

#define EL_IS_ZERO is_empty

#undef ZERO

#define ZERO zero

#undef IS_ZERO

#define IS_ZERO is_zero

#undef FIELD

#define FIELD fold

#undef DEFAULT_IS_ZERO

#define DEFAULT_IS_ZERO 1

#define NO_NEG

#define NO_PULLBACK

#include <isl_pw_templ.c>

#undef UNION

#define UNION isl_union_pw_qpolynomial_fold

#undef PART

#define PART isl_pw_qpolynomial_fold

#undef PARTS

#define PARTS pw_qpolynomial_fold

#define ALIGN_DOMAIN

#define NO_SUB

#include <isl_union_templ.c>

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,

	__isl_take isl_space *dim)

{

	return qpolynomial_fold_alloc(type, dim, 0);

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(

	enum isl_fold type, __isl_take isl_qpolynomial *qp)

{

	isl_qpolynomial_fold *fold;

	if (!qp)

		return NULL;

	fold = qpolynomial_fold_alloc(type, isl_space_copy(qp->dim), 1);

	if (!fold)

		goto error;

	fold->qp[0] = qp;

	fold->n++;

	return fold;

error:

	isl_qpolynomial_fold_free(fold);

	isl_qpolynomial_free(qp);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(

	__isl_keep isl_qpolynomial_fold *fold)

{

	if (!fold)

		return NULL;

	fold->ref++;

	return fold;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(

	__isl_keep isl_qpolynomial_fold *fold)

{

	int i;

	isl_qpolynomial_fold *dup;

	if (!fold)

		return NULL;

	dup = qpolynomial_fold_alloc(fold->type,

					isl_space_copy(fold->dim), fold->n);

	if (!dup)

		return NULL;

	dup->n = fold->n;

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

		dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);

		if (!dup->qp[i])

			goto error;

	}

	return dup;

error:

	isl_qpolynomial_fold_free(dup);

	return NULL;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(

	__isl_take isl_qpolynomial_fold *fold)

{

	if (!fold)

		return NULL;

	if (fold->ref == 1)

		return fold;

	fold->ref--;

	return isl_qpolynomial_fold_dup(fold);

}

void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)

{

	int i;

	if (!fold)

		return;

	if (--fold->ref > 0)

		return;

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

		isl_qpolynomial_free(fold->qp[i]);

	isl_space_free(fold->dim);

	free(fold);

}

int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)

{

	if (!fold)

		return -1;

	return fold->n == 0;

}

__isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(

	__isl_take isl_qpolynomial_fold *fold1,

	__isl_take isl_qpolynomial_fold *fold2)

{

	int i;

	struct isl_qpolynomial_fold *res = NULL;

	if (!fold1 || !fold2)

		goto error;

	isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);

	isl_assert(fold1->dim->ctx, isl_space_is_equal(fold1->dim, fold2->dim),

			goto error);

	if (isl_qpolynomial_fold_is_empty(fold1)) {

		isl_qpolynomial_fold_free(fold1);

		return fold2;

	}

	if (isl_qpolynomial_fold_is_empty(fold2)) {

		isl_qpolynomial_fold_free(fold2);

		return fold1;

	}

	res = qpolynomial_fold_alloc(fold1->type, isl_space_copy(fold1->dim),

					fold1->n + fold2->n);

	if (!res)

		goto error;

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

		res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);

		if (!res->qp[res->n])

			goto error;

		res->n++;

	}

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

		res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);

		if (!res->qp[res->n])

			goto error;

		res->n++;

	}

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return res;

error:

	isl_qpolynomial_fold_free(res);

	isl_qpolynomial_fold_free(fold1);

	isl_qpolynomial_fold_free(fold2);

	return NULL;

}

__isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_fold(

	__isl_take isl_pw_qpolynomial_fold *pw1,

	__isl_take isl_pw_qpolynomial_fold *pw2)

{

	int i, j, n;

	struct isl_pw_qpolynomial_fold *res;

	isl_set *set;

	if (!pw1 || !pw2)

		goto error;

	isl_assert(pw1->dim->ctx, isl_space_is_equal(pw1->dim, pw2->dim), goto error);

	if (isl_pw_qpolynomial_fold_is_zero(pw1)) {

		isl_pw_qpolynomial_fold_free(pw1);

		return pw2;

	}

	if (isl_pw_qpolynomial_fold_is_zero(pw2)) {

		isl_pw_qpolynomial_fold_free(pw2);

		return pw1;

	}

	if (pw1->type != pw2->type)

		isl_die(pw1->dim->ctx, isl_error_invalid,

			"fold types don't match", goto error);

	n = (pw1->n + 1) * (pw2->n + 1);

	res = isl_pw_qpolynomial_fold_alloc_size(isl_space_copy(pw1->dim),

						pw1->type, n);

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

		set = isl_set_copy(pw1->p[i].set);

		for (j = 0; j < pw2->n; ++j) {

			struct isl_set *common;

			isl_qpolynomial_fold *sum;

			set = isl_set_subtract(set,

					isl_set_copy(pw2->p[j].set));

			common = isl_set_intersect(isl_set_copy(pw1->p[i].set),

						isl_set_copy(pw2->p[j].set));

			if (isl_set_plain_is_empty(common)) {

				isl_set_free(common);

				continue;

			}

			sum = isl_qpolynomial_fold_fold_on_domain(common,

			       isl_qpolynomial_fold_copy(pw1->p[i].fold),

			       isl_qpolynomial_fold_copy(pw2->p[j].fold));

			res = isl_pw_qpolynomial_fold_add_piece(res, common, sum);

		}

		res = isl_pw_qpolynomial_fold_add_piece(res, set,

			isl_qpolynomial_fold_copy(pw1->p[i].fold));

	}

	for (j = 0; j < pw2->n; ++j) {

		set = isl_set_copy(pw2->p[j].set);

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

			set = isl_set_subtract(set, isl_set_copy(pw1->p[i].set));

		res = isl_pw_qpolynomial_fold_add_piece(res, set,

				    isl_qpolynomial_fold_copy(pw2->p[j].fold));

	}

	isl_pw_qpolynomial_fold_free(pw1);

	isl_pw_qpolynomial_fold_free(pw2);

	return res;

error: