X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_polynomial.c;h=db081ee830d3dcf555e0d5b9da31f4c8c2ce3767;hb=0b90c0e67236a9f47530297a3fa45dfb2cb06a7a;hp=03d7f1ee58a95ccfeef2aea576462fbced3b9adb;hpb=f22467bba9d54c138441bc60b1cd4b9976bb1804;p=platform%2Fupstream%2Fisl.git

diff --git a/isl_polynomial.c b/isl_polynomial.c
index 03d7f1e..db081ee 100644
--- a/isl_polynomial.c
+++ b/isl_polynomial.c
@@ -9,12 +9,19 @@
  */
 
 #include <stdlib.h>
-#include <isl_seq.h>
+#include <isl_ctx_private.h>
+#include <isl_map_private.h>
+#include <isl_factorization.h>
+#include <isl/lp.h>
+#include <isl/seq.h>
 #include <isl_union_map_private.h>
 #include <isl_polynomial_private.h>
 #include <isl_point_private.h>
 #include <isl_dim_private.h>
-#include <isl_map_private.h>
+#include <isl_div_private.h>
+#include <isl_mat_private.h>
+#include <isl_range.h>
+#include <isl_local_space_private.h>
 
 static unsigned pos(__isl_keep isl_dim *dim, enum isl_dim_type type)
 {
@@ -238,6 +245,20 @@ __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
 	return &cst->up;
 }
 
+__isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
+{
+	struct isl_upoly_cst *cst;
+
+	cst = isl_upoly_cst_alloc(ctx);
+	if (!cst)
+		return NULL;
+
+	isl_int_set_si(cst->n, 1);
+	isl_int_set_si(cst->d, 1);
+
+	return &cst->up;
+}
+
 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
 {
 	struct isl_upoly_cst *cst;
@@ -656,7 +677,50 @@ error:
 	return NULL;
 }
 
-__isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
+__isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
+	__isl_take struct isl_upoly *up, isl_int v)
+{
+	struct isl_upoly_cst *cst;
+
+	up = isl_upoly_cow(up);
+	if (!up)
+		return NULL;
+
+	cst = isl_upoly_as_cst(up);
+
+	isl_int_addmul(cst->n, cst->d, v);
+
+	return up;
+}
+
+__isl_give struct isl_upoly *isl_upoly_add_isl_int(
+	__isl_take struct isl_upoly *up, isl_int v)
+{
+	struct isl_upoly_rec *rec;
+
+	if (!up)
+		return NULL;
+
+	if (isl_upoly_is_cst(up))
+		return isl_upoly_cst_add_isl_int(up, v);
+
+	up = isl_upoly_cow(up);
+	rec = isl_upoly_as_rec(up);
+	if (!rec)
+		goto error;
+
+	rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
+	if (!rec->p[0])
+		goto error;
+
+	return up;
+error:
+	isl_upoly_free(up);
+	return NULL;
+}
+
+__isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
+	__isl_take struct isl_upoly *up, isl_int v)
 {
 	struct isl_upoly_cst *cst;
 
@@ -669,12 +733,13 @@ __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
 
 	cst = isl_upoly_as_cst(up);
 
-	isl_int_neg(cst->n, cst->n);
+	isl_int_mul(cst->n, cst->n, v);
 
 	return up;
 }
 
-__isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
+__isl_give struct isl_upoly *isl_upoly_mul_isl_int(
+	__isl_take struct isl_upoly *up, isl_int v)
 {
 	int i;
 	struct isl_upoly_rec *rec;
@@ -683,7 +748,7 @@ __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
 		return NULL;
 
 	if (isl_upoly_is_cst(up))
-		return isl_upoly_neg_cst(up);
+		return isl_upoly_cst_mul_isl_int(up, v);
 
 	up = isl_upoly_cow(up);
 	rec = isl_upoly_as_rec(up);
@@ -691,7 +756,7 @@ __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
 		goto error;
 
 	for (i = 0; i < rec->n; ++i) {
-		rec->p[i] = isl_upoly_neg(rec->p[i]);
+		rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
 		if (!rec->p[i])
 			goto error;
 	}
@@ -854,6 +919,31 @@ error:
 	return NULL;
 }
 
+__isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
+	unsigned power)
+{
+	struct isl_upoly *res;
+
+	if (!up)
+		return NULL;
+	if (power == 1)
+		return up;
+
+	if (power % 2)
+		res = isl_upoly_copy(up);
+	else
+		res = isl_upoly_one(up->ctx);
+
+	while (power >>= 1) {
+		up = isl_upoly_mul(up, isl_upoly_copy(up));
+		if (power % 2)
+			res = isl_upoly_mul(res, isl_upoly_copy(up));
+	}
+
+	isl_upoly_free(up);
+	return res;
+}
+
 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
 	unsigned n_div, __isl_take struct isl_upoly *up)
 {
@@ -942,7 +1032,7 @@ void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
 	free(qp);
 }
 
-__isl_give struct isl_upoly *isl_upoly_pow(isl_ctx *ctx, int pos, int power)
+__isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
 {
 	int i;
 	struct isl_upoly *up;
@@ -986,7 +1076,7 @@ static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
 
 	isl_assert(up->ctx, rec->n >= 1, goto error);
 
-	base = isl_upoly_pow(up->ctx, r[up->var], 1);
+	base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
 	res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
 
 	for (i = rec->n - 2; i >= 0; --i) {
@@ -1027,19 +1117,6 @@ static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
 	return equal;
 }
 
-static void expand_row(__isl_keep isl_mat *dst, int d,
-	__isl_keep isl_mat *src, int s, int *exp)
-{
-	int i;
-	unsigned c = src->n_col - src->n_row;
-
-	isl_seq_cpy(dst->row[d], src->row[s], c);
-	isl_seq_clr(dst->row[d] + c, dst->n_col - c);
-
-	for (i = 0; i < s; ++i)
-		isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
-}
-
 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
 {
 	int li, lj;
@@ -1067,11 +1144,15 @@ static int div_sort_cmp(const void *p1, const void *p2)
 	return cmp_row(i1->div, i1->row, i2->row);
 }
 
+/* Sort divs and remove duplicates.
+ */
 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
 {
 	int i;
+	int skip;
+	int len;
 	struct isl_div_sort_info *array = NULL;
-	int *pos = NULL;
+	int *pos = NULL, *at = NULL;
 	int *reordering = NULL;
 	unsigned div_pos;
 
@@ -1085,14 +1166,17 @@ static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
 	array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
 				qp->div->n_row);
 	pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
-	reordering = isl_alloc_array(qp->div->ctx, int, qp->div->n_col - 2);
-	if (!array || !pos || !reordering)
+	at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
+	len = qp->div->n_col - 2;
+	reordering = isl_alloc_array(qp->div->ctx, int, len);
+	if (!array || !pos || !at || !reordering)
 		goto error;
 
 	for (i = 0; i < qp->div->n_row; ++i) {
 		array[i].div = qp->div;
 		array[i].row = i;
 		pos[i] = i;
+		at[i] = i;
 	}
 
 	qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
@@ -1101,18 +1185,29 @@ static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
 	for (i = 0; i < div_pos; ++i)
 		reordering[i] = i;
 
-	for (i = 0; i < qp->div->n_row; ++i)
-		reordering[div_pos + array[i].row] = div_pos + i;
-
 	for (i = 0; i < qp->div->n_row; ++i) {
-		int t;
 		if (pos[array[i].row] == i)
 			continue;
-		qp->div = isl_mat_cow(qp->div);
 		qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
-		t = pos[array[i].row];
-		pos[array[i].row] = pos[i];
-		pos[i] = t;
+		pos[at[i]] = pos[array[i].row];
+		at[pos[array[i].row]] = at[i];
+		at[i] = array[i].row;
+		pos[array[i].row] = i;
+	}
+
+	skip = 0;
+	for (i = 0; i < len - div_pos; ++i) {
+		if (i > 0 &&
+		    isl_seq_eq(qp->div->row[i - skip - 1],
+			       qp->div->row[i - skip], qp->div->n_col)) {
+			qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
+			isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
+						 2 + div_pos + i - skip);
+			qp->div = isl_mat_drop_cols(qp->div,
+						    2 + div_pos + i - skip, 1);
+			skip++;
+		}
+		reordering[div_pos + array[i].row] = div_pos + i - skip;
 	}
 
 	qp->upoly = reorder(qp->upoly, reordering);
@@ -1120,12 +1215,14 @@ static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
 	if (!qp->upoly || !qp->div)
 		goto error;
 
+	free(at);
 	free(pos);
 	free(array);
 	free(reordering);
 
 	return qp;
 error:
+	free(at);
 	free(pos);
 	free(array);
 	free(reordering);
@@ -1133,50 +1230,6 @@ error:
 	return NULL;
 }
 
-static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
-	__isl_keep isl_mat *div2, int *exp1, int *exp2)
-{
-	int i, j, k;
-	isl_mat *div = NULL;
-	unsigned d = div1->n_col - div1->n_row;
-
-	div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
-				d + div1->n_row + div2->n_row);
-	if (!div)
-		return NULL;
-
-	for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
-		int cmp;
-
-		expand_row(div, k, div1, i, exp1);
-		expand_row(div, k + 1, div2, j, exp2);
-
-		cmp = cmp_row(div, k, k + 1);
-		if (cmp == 0) {
-			exp1[i++] = k;
-			exp2[j++] = k;
-		} else if (cmp < 0) {
-			exp1[i++] = k;
-		} else {
-			exp2[j++] = k;
-			isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
-		}
-	}
-	for (; i < div1->n_row; ++i, ++k) {
-		expand_row(div, k, div1, i, exp1);
-		exp1[i] = k;
-	}
-	for (; j < div2->n_row; ++j, ++k) {
-		expand_row(div, k, div2, j, exp2);
-		exp2[j] = k;
-	}
-
-	div->n_row = k;
-	div->n_col = d + k;
-
-	return div;
-}
-
 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
 	int *exp, int first)
 {
@@ -1237,7 +1290,7 @@ static __isl_give isl_qpolynomial *with_merged_divs(
 	if (!exp1 || !exp2)
 		goto error;
 
-	div = merge_divs(qp1->div, qp2->div, exp1, exp2);
+	div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
 	if (!div)
 		goto error;
 
@@ -1299,7 +1352,9 @@ __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
 	__isl_take isl_qpolynomial *qp1,
 	__isl_take isl_qpolynomial *qp2)
 {
-	return isl_qpolynomial_add(qp1, qp2);
+	qp1 = isl_qpolynomial_add(qp1, qp2);
+	qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
+	return qp1;
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
@@ -1308,14 +1363,53 @@ __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
 	return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
 }
 
-__isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
+__isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
+	__isl_take isl_qpolynomial *qp, isl_int v)
 {
+	if (isl_int_is_zero(v))
+		return qp;
+
 	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		return NULL;
+
+	qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
+	if (!qp->upoly)
+		goto error;
+
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	return NULL;
+
+}
+
+__isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
+{
+	if (!qp)
+		return NULL;
+
+	return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
+}
 
+__isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
+	__isl_take isl_qpolynomial *qp, isl_int v)
+{
+	if (isl_int_is_one(v))
+		return qp;
+
+	if (qp && isl_int_is_zero(v)) {
+		isl_qpolynomial *zero;
+		zero = isl_qpolynomial_zero(isl_dim_copy(qp->dim));
+		isl_qpolynomial_free(qp);
+		return zero;
+	}
+	
+	qp = isl_qpolynomial_cow(qp);
 	if (!qp)
 		return NULL;
 
-	qp->upoly = isl_upoly_neg(qp->upoly);
+	qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
 	if (!qp->upoly)
 		goto error;
 
@@ -1353,11 +1447,34 @@ error:
 	return NULL;
 }
 
+__isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
+	unsigned power)
+{
+	qp = isl_qpolynomial_cow(qp);
+
+	if (!qp)
+		return NULL;
+
+	qp->upoly = isl_upoly_pow(qp->upoly, power);
+	if (!qp->upoly)
+		goto error;
+
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	return NULL;
+}
+
 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
 {
 	return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
 }
 
+__isl_give isl_qpolynomial *isl_qpolynomial_one(__isl_take isl_dim *dim)
+{
+	return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
+}
+
 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
 {
 	return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
@@ -1513,13 +1630,12 @@ __isl_give isl_vec *isl_qpolynomial_extract_affine(
 	if (!qp)
 		return NULL;
 
-	isl_assert(qp->div->ctx, qp->div->n_row == 0, return NULL);
 	d = isl_dim_total(qp->dim);
-	aff = isl_vec_alloc(qp->div->ctx, 2 + d);
+	aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
 	if (!aff)
 		return NULL;
 
-	isl_seq_clr(aff->el + 1, 1 + d);
+	isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
 	isl_int_set_si(aff->el[0], 1);
 
 	if (isl_upoly_update_affine(qp->upoly, aff) < 0)
@@ -1570,7 +1686,7 @@ void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
 	upoly_update_den(qp->upoly, d);
 }
 
-__isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
+__isl_give isl_qpolynomial *isl_qpolynomial_var_pow(__isl_take isl_dim *dim,
 	int pos, int power)
 {
 	struct isl_ctx *ctx;
@@ -1580,7 +1696,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
 
 	ctx = dim->ctx;
 
-	return isl_qpolynomial_alloc(dim, 0, isl_upoly_pow(ctx, pos, power));
+	return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
@@ -1595,108 +1711,428 @@ __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
 	if (type == isl_dim_set)
 		pos += isl_dim_size(dim, isl_dim_param);
 
-	return isl_qpolynomial_pow(dim, pos, 1);
+	return isl_qpolynomial_var_pow(dim, pos, 1);
 error:
 	isl_dim_free(dim);
 	return NULL;
 }
 
-__isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
-	int power)
+__isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
+	unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
 {
-	struct isl_qpolynomial *qp = NULL;
-	struct isl_upoly_rec *rec;
-	struct isl_upoly_cst *cst;
 	int i;
-	int pos;
+	struct isl_upoly_rec *rec;
+	struct isl_upoly *base, *res;
 
-	if (!div)
+	if (!up)
 		return NULL;
-	isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
 
-	pos = isl_dim_total(div->bmap->dim);
-	rec = isl_upoly_alloc_rec(div->ctx, pos, 1 + power);
-	qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1,
-				   &rec->up);
-	if (!qp)
+	if (isl_upoly_is_cst(up))
+		return up;
+
+	if (up->var < first)
+		return up;
+
+	rec = isl_upoly_as_rec(up);
+	if (!rec)
 		goto error;
 
-	isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
-	isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
+	isl_assert(up->ctx, rec->n >= 1, goto error);
 
-	for (i = 0; i < 1 + power; ++i) {
-		rec->p[i] = isl_upoly_zero(div->ctx);
-		if (!rec->p[i])
-			goto error;
-		rec->n++;
-	}
-	cst = isl_upoly_as_cst(rec->p[power]);
-	isl_int_set_si(cst->n, 1);
+	if (up->var >= first + n)
+		base = isl_upoly_var_pow(up->ctx, up->var, 1);
+	else
+		base = isl_upoly_copy(subs[up->var - first]);
 
-	isl_div_free(div);
+	res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
+	for (i = rec->n - 2; i >= 0; --i) {
+		struct isl_upoly *t;
+		t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
+		res = isl_upoly_mul(res, isl_upoly_copy(base));
+		res = isl_upoly_sum(res, t);
+	}
 
-	return qp;
+	isl_upoly_free(base);
+	isl_upoly_free(up);
+				
+	return res;
 error:
-	isl_qpolynomial_free(qp);
-	isl_div_free(div);
+	isl_upoly_free(up);
 	return NULL;
-}
+}	
 
-__isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
+__isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
+	isl_int denom, unsigned len)
 {
-	return isl_qpolynomial_div_pow(div, 1);
-}
+	int i;
+	struct isl_upoly *up;
 
-__isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
-	const isl_int n, const isl_int d)
-{
-	struct isl_qpolynomial *qp;
-	struct isl_upoly_cst *cst;
+	isl_assert(ctx, len >= 1, return NULL);
 
-	qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
-	if (!qp)
-		return NULL;
+	up = isl_upoly_rat_cst(ctx, f[0], denom);
+	for (i = 0; i < len - 1; ++i) {
+		struct isl_upoly *t;
+		struct isl_upoly *c;
 
-	cst = isl_upoly_as_cst(qp->upoly);
-	isl_int_set(cst->n, n);
-	isl_int_set(cst->d, d);
+		if (isl_int_is_zero(f[1 + i]))
+			continue;
 
-	return qp;
+		c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
+		t = isl_upoly_var_pow(ctx, i, 1);
+		t = isl_upoly_mul(c, t);
+		up = isl_upoly_sum(up, t);
+	}
+
+	return up;
 }
 
-static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
+/* Remove common factor of non-constant terms and denominator.
+ */
+static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
 {
-	struct isl_upoly_rec *rec;
-	int i;
-
-	if (!up)
-		return -1;
-
-	if (isl_upoly_is_cst(up))
-		return 0;
-
-	if (up->var < d)
-		active[up->var] = 1;
+	isl_ctx *ctx = qp->div->ctx;
+	unsigned total = qp->div->n_col - 2;
 
-	rec = isl_upoly_as_rec(up);
-	for (i = 0; i < rec->n; ++i)
-		if (up_set_active(rec->p[i], active, d) < 0)
-			return -1;
+	isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
+	isl_int_gcd(ctx->normalize_gcd,
+		    ctx->normalize_gcd, qp->div->row[div][0]);
+	if (isl_int_is_one(ctx->normalize_gcd))
+		return;
 
-	return 0;
+	isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
+			    ctx->normalize_gcd, total);
+	isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
+			    ctx->normalize_gcd);
+	isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
+			    ctx->normalize_gcd);
 }
 
-static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
+/* Replace the integer division identified by "div" by the polynomial "s".
+ * The integer division is assumed not to appear in the definition
+ * of any other integer divisions.
+ */
+static __isl_give isl_qpolynomial *substitute_div(
+	__isl_take isl_qpolynomial *qp,
+	int div, __isl_take struct isl_upoly *s)
 {
-	int i, j;
-	int d = isl_dim_total(qp->dim);
+	int i;
+	int total;
+	int *reordering;
 
-	if (!qp || !active)
-		return -1;
+	if (!qp || !s)
+		goto error;
 
-	for (i = 0; i < d; ++i)
-		for (j = 0; j < qp->div->n_row; ++j) {
-			if (isl_int_is_zero(qp->div->row[j][2 + i]))
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		goto error;
+
+	total = isl_dim_total(qp->dim);
+	qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
+	if (!qp->upoly)
+		goto error;
+
+	reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
+	if (!reordering)
+		goto error;
+	for (i = 0; i < total + div; ++i)
+		reordering[i] = i;
+	for (i = total + div + 1; i < total + qp->div->n_row; ++i)
+		reordering[i] = i - 1;
+	qp->div = isl_mat_drop_rows(qp->div, div, 1);
+	qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
+	qp->upoly = reorder(qp->upoly, reordering);
+	free(reordering);
+
+	if (!qp->upoly || !qp->div)
+		goto error;
+
+	isl_upoly_free(s);
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	isl_upoly_free(s);
+	return NULL;
+}
+
+/* Replace all integer divisions [e/d] that turn out to not actually be integer
+ * divisions because d is equal to 1 by their definition, i.e., e.
+ */
+static __isl_give isl_qpolynomial *substitute_non_divs(
+	__isl_take isl_qpolynomial *qp)
+{
+	int i, j;
+	int total;
+	struct isl_upoly *s;
+
+	if (!qp)
+		return NULL;
+
+	total = isl_dim_total(qp->dim);
+	for (i = 0; qp && i < qp->div->n_row; ++i) {
+		if (!isl_int_is_one(qp->div->row[i][0]))
+			continue;
+		for (j = i + 1; j < qp->div->n_row; ++j) {
+			if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
+				continue;
+			isl_seq_combine(qp->div->row[j] + 1,
+				qp->div->ctx->one, qp->div->row[j] + 1,
+				qp->div->row[j][2 + total + i],
+				qp->div->row[i] + 1, 1 + total + i);
+			isl_int_set_si(qp->div->row[j][2 + total + i], 0);
+			normalize_div(qp, j);
+		}
+		s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
+					qp->div->row[i][0], qp->div->n_col - 1);
+		qp = substitute_div(qp, i, s);
+		--i;
+	}
+
+	return qp;
+}
+
+/* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
+ * with d the denominator.  When replacing the coefficient e of x by
+ * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
+ * inside the division, so we need to add floor(e/d) * x outside.
+ * That is, we replace q by q' + floor(e/d) * x and we therefore need
+ * to adjust the coefficient of x in each later div that depends on the
+ * current div "div" and also in the affine expression "aff"
+ * (if it too depends on "div").
+ */
+static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
+	__isl_keep isl_vec *aff)
+{
+	int i, j;
+	isl_int v;
+	unsigned total = qp->div->n_col - qp->div->n_row - 2;
+
+	isl_int_init(v);
+	for (i = 0; i < 1 + total + div; ++i) {
+		if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
+		    isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
+			continue;
+		isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
+		isl_int_fdiv_r(qp->div->row[div][1 + i],
+				qp->div->row[div][1 + i], qp->div->row[div][0]);
+		if (!isl_int_is_zero(aff->el[1 + total + div]))
+			isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
+		for (j = div + 1; j < qp->div->n_row; ++j) {
+			if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
+				continue;
+			isl_int_addmul(qp->div->row[j][1 + i],
+					v, qp->div->row[j][2 + total + div]);
+		}
+	}
+	isl_int_clear(v);
+}
+
+/* Check if the last non-zero coefficient is bigger that half of the
+ * denominator.  If so, we will invert the div to further reduce the number
+ * of distinct divs that may appear.
+ * If the last non-zero coefficient is exactly half the denominator,
+ * then we continue looking for earlier coefficients that are bigger
+ * than half the denominator.
+ */
+static int needs_invert(__isl_keep isl_mat *div, int row)
+{
+	int i;
+	int cmp;
+
+	for (i = div->n_col - 1; i >= 1; --i) {
+		if (isl_int_is_zero(div->row[row][i]))
+			continue;
+		isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
+		cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
+		isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
+		if (cmp)
+			return cmp > 0;
+		if (i == 1)
+			return 1;
+	}
+
+	return 0;
+}
+
+/* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
+ * We only invert the coefficients of e (and the coefficient of q in
+ * later divs and in "aff").  After calling this function, the
+ * coefficients of e should be reduced again.
+ */
+static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
+	__isl_keep isl_vec *aff)
+{
+	unsigned total = qp->div->n_col - qp->div->n_row - 2;
+
+	isl_seq_neg(qp->div->row[div] + 1,
+		    qp->div->row[div] + 1, qp->div->n_col - 1);
+	isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
+	isl_int_add(qp->div->row[div][1],
+		    qp->div->row[div][1], qp->div->row[div][0]);
+	if (!isl_int_is_zero(aff->el[1 + total + div]))
+		isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
+	isl_mat_col_mul(qp->div, 2 + total + div,
+			qp->div->ctx->negone, 2 + total + div);
+}
+
+/* Assuming "qp" is a monomial, reduce all its divs to have coefficients
+ * in the interval [0, d-1], with d the denominator and such that the
+ * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
+ *
+ * After the reduction, some divs may have become redundant or identical,
+ * so we call substitute_non_divs and sort_divs.  If these functions
+ * eliminate divs of merge * two or more divs into one, the coefficients
+ * of the enclosing divs may have to be reduced again, so we call
+ * ourselves recursively if the number of divs decreases.
+ */
+static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
+{
+	int i, j;
+	isl_vec *aff = NULL;
+	struct isl_upoly *s;
+	unsigned n_div;
+
+	if (!qp)
+		return NULL;
+
+	aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
+	aff = isl_vec_clr(aff);
+	if (!aff)
+		goto error;
+
+	isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
+
+	for (i = 0; i < qp->div->n_row; ++i) {
+		normalize_div(qp, i);
+		reduce_div(qp, i, aff);
+		if (needs_invert(qp->div, i)) {
+			invert_div(qp, i, aff);
+			reduce_div(qp, i, aff);
+		}
+	}
+
+	s = isl_upoly_from_affine(qp->div->ctx, aff->el,
+				  qp->div->ctx->one, aff->size);
+	qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
+	isl_upoly_free(s);
+	if (!qp->upoly)
+		goto error;
+
+	isl_vec_free(aff);
+
+	n_div = qp->div->n_row;
+	qp = substitute_non_divs(qp);
+	qp = sort_divs(qp);
+	if (qp && qp->div->n_row < n_div)
+		return reduce_divs(qp);
+
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	isl_vec_free(aff);
+	return NULL;
+}
+
+/* Assumes each div only depends on earlier divs.
+ */
+__isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
+	int power)
+{
+	struct isl_qpolynomial *qp = NULL;
+	struct isl_upoly_rec *rec;
+	struct isl_upoly_cst *cst;
+	int i, d;
+	int pos;
+
+	if (!div)
+		return NULL;
+
+	d = div->line - div->bmap->div;
+
+	pos = isl_dim_total(div->bmap->dim) + d;
+	rec = isl_upoly_alloc_rec(div->ctx, pos, 1 + power);
+	qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap),
+				   div->bmap->n_div, &rec->up);
+	if (!qp)
+		goto error;
+
+	for (i = 0; i < div->bmap->n_div; ++i)
+		isl_seq_cpy(qp->div->row[i], div->bmap->div[i], qp->div->n_col);
+
+	for (i = 0; i < 1 + power; ++i) {
+		rec->p[i] = isl_upoly_zero(div->ctx);
+		if (!rec->p[i])
+			goto error;
+		rec->n++;
+	}
+	cst = isl_upoly_as_cst(rec->p[power]);
+	isl_int_set_si(cst->n, 1);
+
+	isl_div_free(div);
+
+	qp = reduce_divs(qp);
+
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	isl_div_free(div);
+	return NULL;
+}
+
+__isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
+{
+	return isl_qpolynomial_div_pow(div, 1);
+}
+
+__isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
+	const isl_int n, const isl_int d)
+{
+	struct isl_qpolynomial *qp;
+	struct isl_upoly_cst *cst;
+
+	qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
+	if (!qp)
+		return NULL;
+
+	cst = isl_upoly_as_cst(qp->upoly);
+	isl_int_set(cst->n, n);
+	isl_int_set(cst->d, d);
+
+	return qp;
+}
+
+static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
+{
+	struct isl_upoly_rec *rec;
+	int i;
+
+	if (!up)
+		return -1;
+
+	if (isl_upoly_is_cst(up))
+		return 0;
+
+	if (up->var < d)
+		active[up->var] = 1;
+
+	rec = isl_upoly_as_rec(up);
+	for (i = 0; i < rec->n; ++i)
+		if (up_set_active(rec->p[i], active, d) < 0)
+			return -1;
+
+	return 0;
+}
+
+static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
+{
+	int i, j;
+	int d = isl_dim_total(qp->dim);
+
+	if (!qp || !active)
+		return -1;
+
+	for (i = 0; i < d; ++i)
+		for (j = 0; j < qp->div->n_row; ++j) {
+			if (isl_int_is_zero(qp->div->row[j][2 + i]))
 				continue;
 			active[i] = 1;
 			break;
@@ -1776,6 +2212,22 @@ error:
 	return NULL;
 }
 
+__isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
+	__isl_take isl_qpolynomial *qp,
+	enum isl_dim_type type, unsigned pos, const char *s)
+{
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		return NULL;
+	qp->dim = isl_dim_set_name(qp->dim, type, pos, s);
+	if (!qp->dim)
+		goto error;
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	return NULL;
+}
+
 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
 	__isl_take isl_qpolynomial *qp,
 	enum isl_dim_type type, unsigned first, unsigned n)
@@ -1815,112 +2267,44 @@ error:
 	return NULL;
 }
 
-__isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
-	unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
+__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
+	__isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
 {
-	int i;
-	struct isl_upoly_rec *rec;
-	struct isl_upoly *base, *res;
-
-	if (!up)
-		return NULL;
-
-	if (isl_upoly_is_cst(up))
-		return up;
+	int i, j, k;
+	isl_int denom;
+	unsigned total;
+	unsigned n_div;
+	struct isl_upoly *up;
 
-	if (up->var < first)
-		return up;
+	if (!eq)
+		goto error;
+	if (eq->n_eq == 0) {
+		isl_basic_set_free(eq);
+		return qp;
+	}
 
-	rec = isl_upoly_as_rec(up);
-	if (!rec)
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		goto error;
+	qp->div = isl_mat_cow(qp->div);
+	if (!qp->div)
 		goto error;
 
-	isl_assert(up->ctx, rec->n >= 1, goto error);
+	total = 1 + isl_dim_total(eq->dim);
+	n_div = eq->n_div;
+	isl_int_init(denom);
+	for (i = 0; i < eq->n_eq; ++i) {
+		j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
+		if (j < 0 || j == 0 || j >= total)
+			continue;
 
-	if (up->var >= first + n)
-		base = isl_upoly_pow(up->ctx, up->var, 1);
-	else
-		base = isl_upoly_copy(subs[up->var - first]);
-
-	res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
-	for (i = rec->n - 2; i >= 0; --i) {
-		struct isl_upoly *t;
-		t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
-		res = isl_upoly_mul(res, isl_upoly_copy(base));
-		res = isl_upoly_sum(res, t);
-	}
-
-	isl_upoly_free(base);
-	isl_upoly_free(up);
-				
-	return res;
-error:
-	isl_upoly_free(up);
-	return NULL;
-}	
-
-__isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
-	isl_int denom, unsigned len)
-{
-	int i;
-	struct isl_upoly *up;
-
-	isl_assert(ctx, len >= 1, return NULL);
-
-	up = isl_upoly_rat_cst(ctx, f[0], denom);
-	for (i = 0; i < len - 1; ++i) {
-		struct isl_upoly *t;
-		struct isl_upoly *c;
-
-		if (isl_int_is_zero(f[1 + i]))
-			continue;
-
-		c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
-		t = isl_upoly_pow(ctx, i, 1);
-		t = isl_upoly_mul(c, t);
-		up = isl_upoly_sum(up, t);
-	}
-
-	return up;
-}
-
-__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
-	__isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
-{
-	int i, j, k;
-	isl_int denom;
-	unsigned total;
-	struct isl_upoly *up;
-
-	if (!eq)
-		goto error;
-	if (eq->n_eq == 0) {
-		isl_basic_set_free(eq);
-		return qp;
-	}
-
-	qp = isl_qpolynomial_cow(qp);
-	if (!qp)
-		goto error;
-	qp->div = isl_mat_cow(qp->div);
-	if (!qp->div)
-		goto error;
-
-	total = 1 + isl_dim_total(eq->dim);
-	isl_int_init(denom);
-	for (i = 0; i < eq->n_eq; ++i) {
-		j = isl_seq_last_non_zero(eq->eq[i], total);
-		if (j < 0 || j == 0)
-			continue;
-
-		for (k = 0; k < qp->div->n_row; ++k) {
-			if (isl_int_is_zero(qp->div->row[k][1 + j]))
-				continue;
-			isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
-					&qp->div->row[k][0]);
-			isl_seq_normalize(qp->div->ctx,
-					  qp->div->row[k], 1 + total);
-		}
+		for (k = 0; k < qp->div->n_row; ++k) {
+			if (isl_int_is_zero(qp->div->row[k][1 + j]))
+				continue;
+			isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
+					&qp->div->row[k][0]);
+			normalize_div(qp, k);
+		}
 
 		if (isl_int_is_pos(eq->eq[i][j]))
 			isl_seq_neg(eq->eq[i], eq->eq[i], total);
@@ -1939,6 +2323,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
 
 	isl_basic_set_free(eq);
 
+	qp = substitute_non_divs(qp);
 	qp = sort_divs(qp);
 
 	return qp;
@@ -1948,6 +2333,62 @@ error:
 	return NULL;
 }
 
+static __isl_give isl_basic_set *add_div_constraints(
+	__isl_take isl_basic_set *bset, __isl_take isl_mat *div)
+{
+	int i;
+	unsigned total;
+
+	if (!bset || !div)
+		goto error;
+
+	bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
+	if (!bset)
+		goto error;
+	total = isl_basic_set_total_dim(bset);
+	for (i = 0; i < div->n_row; ++i)
+		if (isl_basic_set_add_div_constraints_var(bset,
+				    total - div->n_row + i, div->row[i]) < 0)
+			goto error;
+
+	isl_mat_free(div);
+	return bset;
+error:
+	isl_mat_free(div);
+	isl_basic_set_free(bset);
+	return NULL;
+}
+
+/* Look for equalities among the variables shared by context and qp
+ * and the integer divisions of qp, if any.
+ * The equalities are then used to eliminate variables and/or integer
+ * divisions from qp.
+ */
+__isl_give isl_qpolynomial *isl_qpolynomial_gist(
+	__isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
+{
+	isl_basic_set *aff;
+
+	if (!qp)
+		goto error;
+	if (qp->div->n_row > 0) {
+		isl_basic_set *bset;
+		context = isl_set_add_dims(context, isl_dim_set,
+					    qp->div->n_row);
+		bset = isl_basic_set_universe(isl_set_get_dim(context));
+		bset = add_div_constraints(bset, isl_mat_copy(qp->div));
+		context = isl_set_intersect(context,
+					    isl_set_from_basic_set(bset));
+	}
+
+	aff = isl_set_affine_hull(context);
+	return isl_qpolynomial_substitute_equalities(qp, aff);
+error:
+	isl_qpolynomial_free(qp);
+	isl_set_free(context);
+	return NULL;
+}
+
 #undef PW
 #define PW isl_pw_qpolynomial
 #undef EL
@@ -1976,7 +2417,7 @@ int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
 	if (pwqp->n != -1)
 		return 0;
 
-	if (!isl_set_fast_is_universe(pwqp->p[0].set))
+	if (!isl_set_plain_is_universe(pwqp->p[0].set))
 		return 0;
 
 	return isl_qpolynomial_is_one(pwqp->p[0].qp);
@@ -2025,7 +2466,7 @@ __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
 			struct isl_qpolynomial *prod;
 			common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
 						isl_set_copy(pwqp2->p[j].set));
-			if (isl_set_fast_is_empty(common)) {
+			if (isl_set_plain_is_empty(common)) {
 				isl_set_free(common);
 				continue;
 			}
@@ -2504,32 +2945,6 @@ error:
 	return NULL;
 }
 
-__isl_give isl_basic_set *add_div_constraints(__isl_take isl_basic_set *bset,
-	__isl_take isl_mat *div)
-{
-	int i;
-	unsigned total;
-
-	if (!bset || !div)
-		goto error;
-
-	bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
-	if (!bset)
-		goto error;
-	total = isl_basic_set_total_dim(bset);
-	for (i = 0; i < div->n_row; ++i)
-		if (isl_basic_set_add_div_constraints_var(bset,
-				    total - div->n_row + i, div->row[i]) < 0)
-			goto error;
-
-	isl_mat_free(div);
-	return bset;
-error:
-	isl_mat_free(div);
-	isl_basic_set_free(bset);
-	return NULL;
-}
-
 /* Extend "bset" with extra set dimensions for each integer division
  * in "qp" and then call "fn" with the extended bset and the polynomial
  * that results from replacing each of the integer divisions by the
@@ -2712,7 +3127,7 @@ __isl_give struct isl_upoly *isl_upoly_homogenize(
 	if (isl_upoly_is_cst(up) || up->var < first) {
 		struct isl_upoly *hom;
 
-		hom = isl_upoly_pow(up->ctx, first, target - deg);
+		hom = isl_upoly_var_pow(up->ctx, first, target - deg);
 		if (!hom)
 			goto error;
 		rec = isl_upoly_as_rec(hom);
@@ -3030,7 +3445,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
 		if (!term->pow[i])
 			continue;
 		up = isl_upoly_mul(up,
-			isl_upoly_pow(term->dim->ctx, i, term->pow[i]));
+			isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
 	}
 
 	qp = isl_qpolynomial_alloc(isl_dim_copy(term->dim), term->div->n_row, up);
@@ -3214,6 +3629,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp
 	__isl_take isl_morph *morph)
 {
 	int i;
+	int n_sub;
 	isl_ctx *ctx;
 	struct isl_upoly *up;
 	unsigned n_div;
@@ -3227,17 +3643,24 @@ __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp
 	ctx = qp->dim->ctx;
 	isl_assert(ctx, isl_dim_equal(qp->dim, morph->dom->dim), goto error);
 
-	subs = isl_calloc_array(ctx, struct isl_upoly *, morph->inv->n_row - 1);
+	n_sub = morph->inv->n_row - 1;
+	if (morph->inv->n_row != morph->inv->n_col)
+		n_sub += qp->div->n_row;
+	subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
 	if (!subs)
 		goto error;
 
 	for (i = 0; 1 + i < morph->inv->n_row; ++i)
 		subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
 					morph->inv->row[0][0], morph->inv->n_col);
+	if (morph->inv->n_row != morph->inv->n_col)
+		for (i = 0; i < qp->div->n_row; ++i)
+			subs[morph->inv->n_row - 1 + i] =
+			    isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
 
-	qp->upoly = isl_upoly_subs(qp->upoly, 0, morph->inv->n_row - 1, subs);
+	qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
 
-	for (i = 0; 1 + i < morph->inv->n_row; ++i)
+	for (i = 0; i < n_sub; ++i)
 		isl_upoly_free(subs[i]);
 	free(subs);
 
@@ -3332,3 +3755,782 @@ __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
 {
 	return match_bin_op(upwqp1, upwqp2, &mul_entry);
 }
+
+/* Reorder the columns of the given div definitions according to the
+ * given reordering.
+ */
+static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
+	__isl_take isl_reordering *r)
+{
+	int i, j;
+	isl_mat *mat;
+	int extra;
+
+	if (!div || !r)
+		goto error;
+
+	extra = isl_dim_total(r->dim) + div->n_row - r->len;
+	mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
+	if (!mat)
+		goto error;
+
+	for (i = 0; i < div->n_row; ++i) {
+		isl_seq_cpy(mat->row[i], div->row[i], 2);
+		isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
+		for (j = 0; j < r->len; ++j)
+			isl_int_set(mat->row[i][2 + r->pos[j]],
+				    div->row[i][2 + j]);
+	}
+
+	isl_reordering_free(r);
+	isl_mat_free(div);
+	return mat;
+error:
+	isl_reordering_free(r);
+	isl_mat_free(div);
+	return NULL;
+}
+
+/* Reorder the dimension of "qp" according to the given reordering.
+ */
+__isl_give isl_qpolynomial *isl_qpolynomial_realign(
+	__isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
+{
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		goto error;
+
+	r = isl_reordering_extend(r, qp->div->n_row);
+	if (!r)
+		goto error;
+
+	qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
+	if (!qp->div)
+		goto error;
+
+	qp->upoly = reorder(qp->upoly, r->pos);
+	if (!qp->upoly)
+		goto error;
+
+	qp = isl_qpolynomial_reset_dim(qp, isl_dim_copy(r->dim));
+
+	isl_reordering_free(r);
+	return qp;
+error:
+	isl_qpolynomial_free(qp);
+	isl_reordering_free(r);
+	return NULL;
+}
+
+__isl_give isl_qpolynomial *isl_qpolynomial_align_params(
+	__isl_take isl_qpolynomial *qp, __isl_take isl_dim *model)
+{
+	if (!qp || !model)
+		goto error;
+
+	if (!isl_dim_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
+		isl_reordering *exp;
+
+		model = isl_dim_drop(model, isl_dim_in,
+					0, isl_dim_size(model, isl_dim_in));
+		model = isl_dim_drop(model, isl_dim_out,
+					0, isl_dim_size(model, isl_dim_out));
+		exp = isl_parameter_alignment_reordering(qp->dim, model);
+		exp = isl_reordering_extend_dim(exp,
+						isl_qpolynomial_get_dim(qp));
+		qp = isl_qpolynomial_realign(qp, exp);
+	}
+
+	isl_dim_free(model);
+	return qp;
+error:
+	isl_dim_free(model);
+	isl_qpolynomial_free(qp);
+	return NULL;
+}
+
+struct isl_split_periods_data {
+	int max_periods;
+	isl_pw_qpolynomial *res;
+};
+
+/* Create a slice where the integer division "div" has the fixed value "v".
+ * In particular, if "div" refers to floor(f/m), then create a slice
+ *
+ *	m v <= f <= m v + (m - 1)
+ *
+ * or
+ *
+ *	f - m v >= 0
+ *	-f + m v + (m - 1) >= 0
+ */
+static __isl_give isl_set *set_div_slice(__isl_take isl_dim *dim,
+	__isl_keep isl_qpolynomial *qp, int div, isl_int v)
+{
+	int total;
+	isl_basic_set *bset = NULL;
+	int k;
+
+	if (!dim || !qp)
+		goto error;
+
+	total = isl_dim_total(dim);
+	bset = isl_basic_set_alloc_dim(isl_dim_copy(dim), 0, 0, 2);
+
+	k = isl_basic_set_alloc_inequality(bset);
+	if (k < 0)
+		goto error;
+	isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
+	isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
+
+	k = isl_basic_set_alloc_inequality(bset);
+	if (k < 0)
+		goto error;
+	isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
+	isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
+	isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
+	isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
+
+	isl_dim_free(dim);
+	return isl_set_from_basic_set(bset);
+error:
+	isl_basic_set_free(bset);
+	isl_dim_free(dim);
+	return NULL;
+}
+
+static int split_periods(__isl_take isl_set *set,
+	__isl_take isl_qpolynomial *qp, void *user);
+
+/* Create a slice of the domain "set" such that integer division "div"
+ * has the fixed value "v" and add the results to data->res,
+ * replacing the integer division by "v" in "qp".
+ */
+static int set_div(__isl_take isl_set *set,
+	__isl_take isl_qpolynomial *qp, int div, isl_int v,
+	struct isl_split_periods_data *data)
+{
+	int i;
+	int total;
+	isl_set *slice;
+	struct isl_upoly *cst;
+
+	slice = set_div_slice(isl_set_get_dim(set), qp, div, v);
+	set = isl_set_intersect(set, slice);
+
+	if (!qp)
+		goto error;
+
+	total = isl_dim_total(qp->dim);
+
+	for (i = div + 1; i < qp->div->n_row; ++i) {
+		if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
+			continue;
+		isl_int_addmul(qp->div->row[i][1],
+				qp->div->row[i][2 + total + div], v);
+		isl_int_set_si(qp->div->row[i][2 + total + div], 0);
+	}
+
+	cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
+	qp = substitute_div(qp, div, cst);
+
+	return split_periods(set, qp, data);
+error:
+	isl_set_free(set);
+	isl_qpolynomial_free(qp);
+	return -1;
+}
+
+/* Split the domain "set" such that integer division "div"
+ * has a fixed value (ranging from "min" to "max") on each slice
+ * and add the results to data->res.
+ */
+static int split_div(__isl_take isl_set *set,
+	__isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
+	struct isl_split_periods_data *data)
+{
+	for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
+		isl_set *set_i = isl_set_copy(set);
+		isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
+
+		if (set_div(set_i, qp_i, div, min, data) < 0)
+			goto error;
+	}
+	isl_set_free(set);
+	isl_qpolynomial_free(qp);
+	return 0;
+error:
+	isl_set_free(set);
+	isl_qpolynomial_free(qp);
+	return -1;
+}
+
+/* If "qp" refers to any integer division
+ * that can only attain "max_periods" distinct values on "set"
+ * then split the domain along those distinct values.
+ * Add the results (or the original if no splitting occurs)
+ * to data->res.
+ */
+static int split_periods(__isl_take isl_set *set,
+	__isl_take isl_qpolynomial *qp, void *user)
+{
+	int i;
+	isl_pw_qpolynomial *pwqp;
+	struct isl_split_periods_data *data;
+	isl_int min, max;
+	int total;
+	int r = 0;
+
+	data = (struct isl_split_periods_data *)user;
+
+	if (!set || !qp)
+		goto error;
+
+	if (qp->div->n_row == 0) {
+		pwqp = isl_pw_qpolynomial_alloc(set, qp);
+		data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
+		return 0;
+	}
+
+	isl_int_init(min);
+	isl_int_init(max);
+	total = isl_dim_total(qp->dim);
+	for (i = 0; i < qp->div->n_row; ++i) {
+		enum isl_lp_result lp_res;
+
+		if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
+						qp->div->n_row) != -1)
+			continue;
+
+		lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
+					  set->ctx->one, &min, NULL, NULL);
+		if (lp_res == isl_lp_error)
+			goto error2;
+		if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
+			continue;
+		isl_int_fdiv_q(min, min, qp->div->row[i][0]);
+
+		lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
+					  set->ctx->one, &max, NULL, NULL);
+		if (lp_res == isl_lp_error)
+			goto error2;
+		if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
+			continue;
+		isl_int_fdiv_q(max, max, qp->div->row[i][0]);
+
+		isl_int_sub(max, max, min);
+		if (isl_int_cmp_si(max, data->max_periods) < 0) {
+			isl_int_add(max, max, min);
+			break;
+		}
+	}
+
+	if (i < qp->div->n_row) {
+		r = split_div(set, qp, i, min, max, data);
+	} else {
+		pwqp = isl_pw_qpolynomial_alloc(set, qp);
+		data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
+	}
+
+	isl_int_clear(max);
+	isl_int_clear(min);
+
+	return r;
+error2:
+	isl_int_clear(max);
+	isl_int_clear(min);
+error:
+	isl_set_free(set);
+	isl_qpolynomial_free(qp);
+	return -1;
+}
+
+/* If any quasi-polynomial in pwqp refers to any integer division
+ * that can only attain "max_periods" distinct values on its domain
+ * then split the domain along those distinct values.
+ */
+__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
+	__isl_take isl_pw_qpolynomial *pwqp, int max_periods)
+{
+	struct isl_split_periods_data data;
+
+	data.max_periods = max_periods;
+	data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_dim(pwqp));
+
+	if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
+		goto error;
+
+	isl_pw_qpolynomial_free(pwqp);
+
+	return data.res;
+error:
+	isl_pw_qpolynomial_free(data.res);
+	isl_pw_qpolynomial_free(pwqp);
+	return NULL;
+}
+
+/* Construct a piecewise quasipolynomial that is constant on the given
+ * domain.  In particular, it is
+ *	0	if cst == 0
+ *	1	if cst == 1
+ *  infinity	if cst == -1
+ */
+static __isl_give isl_pw_qpolynomial *constant_on_domain(
+	__isl_take isl_basic_set *bset, int cst)
+{
+	isl_dim *dim;
+	isl_qpolynomial *qp;
+
+	if (!bset)
+		return NULL;
+
+	bset = isl_basic_map_domain(isl_basic_map_from_range(bset));
+	dim = isl_basic_set_get_dim(bset);
+	if (cst < 0)
+		qp = isl_qpolynomial_infty(dim);
+	else if (cst == 0)
+		qp = isl_qpolynomial_zero(dim);
+	else
+		qp = isl_qpolynomial_one(dim);
+	return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
+}
+
+/* Factor bset, call fn on each of the factors and return the product.
+ *
+ * If no factors can be found, simply call fn on the input.
+ * Otherwise, construct the factors based on the factorizer,
+ * call fn on each factor and compute the product.
+ */
+static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
+	__isl_take isl_basic_set *bset,
+	__isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
+{
+	int i, n;
+	isl_dim *dim;
+	isl_set *set;
+	isl_factorizer *f;
+	isl_qpolynomial *qp;
+	isl_pw_qpolynomial *pwqp;
+	unsigned nparam;
+	unsigned nvar;
+
+	f = isl_basic_set_factorizer(bset);
+	if (!f)
+		goto error;
+	if (f->n_group == 0) {
+		isl_factorizer_free(f);
+		return fn(bset);
+	}
+
+	nparam = isl_basic_set_dim(bset, isl_dim_param);
+	nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+	dim = isl_basic_set_get_dim(bset);
+	dim = isl_dim_domain(dim);
+	set = isl_set_universe(isl_dim_copy(dim));
+	qp = isl_qpolynomial_one(dim);
+	pwqp = isl_pw_qpolynomial_alloc(set, qp);
+
+	bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
+
+	for (i = 0, n = 0; i < f->n_group; ++i) {
+		isl_basic_set *bset_i;
+		isl_pw_qpolynomial *pwqp_i;
+
+		bset_i = isl_basic_set_copy(bset);
+		bset_i = isl_basic_set_drop_constraints_involving(bset_i,
+			    nparam + n + f->len[i], nvar - n - f->len[i]);
+		bset_i = isl_basic_set_drop_constraints_involving(bset_i,
+			    nparam, n);
+		bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
+			    n + f->len[i], nvar - n - f->len[i]);
+		bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
+
+		pwqp_i = fn(bset_i);
+		pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
+
+		n += f->len[i];
+	}
+
+	isl_basic_set_free(bset);
+	isl_factorizer_free(f);
+
+	return pwqp;
+error:
+	isl_basic_set_free(bset);
+	return NULL;
+}
+
+/* Factor bset, call fn on each of the factors and return the product.
+ * The function is assumed to evaluate to zero on empty domains,
+ * to one on zero-dimensional domains and to infinity on unbounded domains
+ * and will not be called explicitly on zero-dimensional or unbounded domains.
+ *
+ * We first check for some special cases and remove all equalities.
+ * Then we hand over control to compressed_multiplicative_call.
+ */
+__isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
+	__isl_take isl_basic_set *bset,
+	__isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
+{
+	int bounded;
+	isl_morph *morph;
+	isl_pw_qpolynomial *pwqp;
+	unsigned orig_nvar, final_nvar;
+
+	if (!bset)
+		return NULL;
+
+	if (isl_basic_set_plain_is_empty(bset))
+		return constant_on_domain(bset, 0);
+
+	orig_nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+	if (orig_nvar == 0)
+		return constant_on_domain(bset, 1);
+
+	bounded = isl_basic_set_is_bounded(bset);
+	if (bounded < 0)
+		goto error;
+	if (!bounded)
+		return constant_on_domain(bset, -1);
+
+	if (bset->n_eq == 0)
+		return compressed_multiplicative_call(bset, fn);
+
+	morph = isl_basic_set_full_compression(bset);
+	bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
+
+	final_nvar = isl_basic_set_dim(bset, isl_dim_set);
+
+	pwqp = compressed_multiplicative_call(bset, fn);
+
+	morph = isl_morph_remove_dom_dims(morph, isl_dim_set, 0, orig_nvar);
+	morph = isl_morph_remove_ran_dims(morph, isl_dim_set, 0, final_nvar);
+	morph = isl_morph_inverse(morph);
+
+	pwqp = isl_pw_qpolynomial_morph(pwqp, morph);
+
+	return pwqp;
+error:
+	isl_basic_set_free(bset);
+	return NULL;
+}
+
+/* Drop all floors in "qp", turning each integer division [a/m] into
+ * a rational division a/m.  If "down" is set, then the integer division
+ * is replaces by (a-(m-1))/m instead.
+ */
+static __isl_give isl_qpolynomial *qp_drop_floors(
+	__isl_take isl_qpolynomial *qp, int down)
+{
+	int i;
+	struct isl_upoly *s;
+
+	if (!qp)
+		return NULL;
+	if (qp->div->n_row == 0)
+		return qp;
+
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		return NULL;
+
+	for (i = qp->div->n_row - 1; i >= 0; --i) {
+		if (down) {
+			isl_int_sub(qp->div->row[i][1],
+				    qp->div->row[i][1], qp->div->row[i][0]);
+			isl_int_add_ui(qp->div->row[i][1],
+				       qp->div->row[i][1], 1);
+		}
+		s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
+					qp->div->row[i][0], qp->div->n_col - 1);
+		qp = substitute_div(qp, i, s);
+		if (!qp)
+			return NULL;
+	}
+
+	return qp;
+}
+
+/* Drop all floors in "pwqp", turning each integer division [a/m] into
+ * a rational division a/m.
+ */
+static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
+	__isl_take isl_pw_qpolynomial *pwqp)
+{
+	int i;
+
+	if (!pwqp)
+		return NULL;
+
+	if (isl_pw_qpolynomial_is_zero(pwqp))
+		return pwqp;
+
+	pwqp = isl_pw_qpolynomial_cow(pwqp);
+	if (!pwqp)
+		return NULL;
+
+	for (i = 0; i < pwqp->n; ++i) {
+		pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
+		if (!pwqp->p[i].qp)
+			goto error;
+	}
+
+	return pwqp;
+error:
+	isl_pw_qpolynomial_free(pwqp);
+	return NULL;
+}
+
+/* Adjust all the integer divisions in "qp" such that they are at least
+ * one over the given orthant (identified by "signs").  This ensures
+ * that they will still be non-negative even after subtracting (m-1)/m.
+ *
+ * In particular, f is replaced by f' + v, changing f = [a/m]
+ * to f' = [(a - m v)/m].
+ * If the constant term k in a is smaller than m,
+ * the constant term of v is set to floor(k/m) - 1.
+ * For any other term, if the coefficient c and the variable x have
+ * the same sign, then no changes are needed.
+ * Otherwise, if the variable is positive (and c is negative),
+ * then the coefficient of x in v is set to floor(c/m).
+ * If the variable is negative (and c is positive),
+ * then the coefficient of x in v is set to ceil(c/m).
+ */
+static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
+	int *signs)
+{
+	int i, j;
+	int total;
+	isl_vec *v = NULL;
+	struct isl_upoly *s;
+
+	qp = isl_qpolynomial_cow(qp);
+	if (!qp)
+		return NULL;
+	qp->div = isl_mat_cow(qp->div);
+	if (!qp->div)
+		goto error;
+
+	total = isl_dim_total(qp->dim);
+	v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
+
+	for (i = 0; i < qp->div->n_row; ++i) {
+		isl_int *row = qp->div->row[i];
+		v = isl_vec_clr(v);
+		if (!v)
+			goto error;
+		if (isl_int_lt(row[1], row[0])) {
+			isl_int_fdiv_q(v->el[0], row[1], row[0]);
+			isl_int_sub_ui(v->el[0], v->el[0], 1);
+			isl_int_submul(row[1], row[0], v->el[0]);
+		}
+		for (j = 0; j < total; ++j) {
+			if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
+				continue;
+			if (signs[j] < 0)
+				isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
+			else
+				isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
+			isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
+		}
+		for (j = 0; j < i; ++j) {
+			if (isl_int_sgn(row[2 + total + j]) >= 0)
+				continue;
+			isl_int_fdiv_q(v->el[1 + total + j],
+					row[2 + total + j], row[0]);
+			isl_int_submul(row[2 + total + j],
+					row[0], v->el[1 + total + j]);
+		}
+		for (j = i + 1; j < qp->div->n_row; ++j) {
+			if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
+				continue;
+			isl_seq_combine(qp->div->row[j] + 1,
+				qp->div->ctx->one, qp->div->row[j] + 1,
+				qp->div->row[j][2 + total + i], v->el, v->size);
+		}
+		isl_int_set_si(v->el[1 + total + i], 1);
+		s = isl_upoly_from_affine(qp->dim->ctx, v->el,
+					qp->div->ctx->one, v->size);
+		qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
+		isl_upoly_free(s);
+		if (!qp->upoly)
+			goto error;
+	}
+
+	isl_vec_free(v);
+	return qp;
+error:
+	isl_vec_free(v);
+	isl_qpolynomial_free(qp);
+	return NULL;
+}
+
+struct isl_to_poly_data {
+	int sign;
+	isl_pw_qpolynomial *res;
+	isl_qpolynomial *qp;
+};
+
+/* Appoximate data->qp by a polynomial on the orthant identified by "signs".
+ * We first make all integer divisions positive and then split the
+ * quasipolynomials into terms with sign data->sign (the direction
+ * of the requested approximation) and terms with the opposite sign.
+ * In the first set of terms, each integer division [a/m] is
+ * overapproximated by a/m, while in the second it is underapproximated
+ * by (a-(m-1))/m.
+ */
+static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
+	void *user)
+{
+	struct isl_to_poly_data *data = user;
+	isl_pw_qpolynomial *t;
+	isl_qpolynomial *qp, *up, *down;
+
+	qp = isl_qpolynomial_copy(data->qp);
+	qp = make_divs_pos(qp, signs);
+
+	up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
+	up = qp_drop_floors(up, 0);
+	down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
+	down = qp_drop_floors(down, 1);
+
+	isl_qpolynomial_free(qp);
+	qp = isl_qpolynomial_add(up, down);
+
+	t = isl_pw_qpolynomial_alloc(orthant, qp);
+	data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
+
+	return 0;
+}
+
+/* Approximate each quasipolynomial by a polynomial.  If "sign" is positive,
+ * the polynomial will be an overapproximation.  If "sign" is negative,
+ * it will be an underapproximation.  If "sign" is zero, the approximation
+ * will lie somewhere in between.
+ *
+ * In particular, is sign == 0, we simply drop the floors, turning
+ * the integer divisions into rational divisions.
+ * Otherwise, we split the domains into orthants, make all integer divisions
+ * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
+ * depending on the requested sign and the sign of the term in which
+ * the integer division appears.
+ */
+__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
+	__isl_take isl_pw_qpolynomial *pwqp, int sign)
+{
+	int i;
+	struct isl_to_poly_data data;
+
+	if (sign == 0)
+		return pwqp_drop_floors(pwqp);
+
+	if (!pwqp)
+		return NULL;
+
+	data.sign = sign;
+	data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_dim(pwqp));
+
+	for (i = 0; i < pwqp->n; ++i) {
+		if (pwqp->p[i].qp->div->n_row == 0) {
+			isl_pw_qpolynomial *t;
+			t = isl_pw_qpolynomial_alloc(
+					isl_set_copy(pwqp->p[i].set),
+					isl_qpolynomial_copy(pwqp->p[i].qp));
+			data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
+			continue;
+		}
+		data.qp = pwqp->p[i].qp;
+		if (isl_set_foreach_orthant(pwqp->p[i].set,
+					&to_polynomial_on_orthant, &data) < 0)
+			goto error;
+	}
+
+	isl_pw_qpolynomial_free(pwqp);
+
+	return data.res;
+error:
+	isl_pw_qpolynomial_free(pwqp);
+	isl_pw_qpolynomial_free(data.res);
+	return NULL;
+}
+
+static int poly_entry(void **entry, void *user)
+{
+	int *sign = user;
+	isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
+
+	*pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
+
+	return *pwqp ? 0 : -1;
+}
+
+__isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
+	__isl_take isl_union_pw_qpolynomial *upwqp, int sign)
+{
+	upwqp = isl_union_pw_qpolynomial_cow(upwqp);
+	if (!upwqp)
+		return NULL;
+
+	if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
+				   &poly_entry, &sign) < 0)
+		goto error;
+
+	return upwqp;
+error:
+	isl_union_pw_qpolynomial_free(upwqp);
+	return NULL;
+}
+
+__isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
+	__isl_take isl_qpolynomial *qp)
+{
+	int i, k;
+	isl_dim *dim;
+	isl_vec *aff = NULL;
+	isl_basic_map *bmap = NULL;
+	unsigned pos;
+	unsigned n_div;
+
+	if (!qp)
+		return NULL;
+	if (!isl_upoly_is_affine(qp->upoly))
+		isl_die(qp->dim->ctx, isl_error_invalid,
+			"input quasi-polynomial not affine", goto error);
+	aff = isl_qpolynomial_extract_affine(qp);
+	if (!aff)
+		goto error;
+	dim = isl_qpolynomial_get_dim(qp);
+	dim = isl_dim_from_domain(dim);
+	pos = 1 + isl_dim_offset(dim, isl_dim_out);
+	dim = isl_dim_add(dim, isl_dim_out, 1);
+	n_div = qp->div->n_row;
+	bmap = isl_basic_map_alloc_dim(dim, n_div, 1, 2 * n_div);
+
+	for (i = 0; i < n_div; ++i) {
+		k = isl_basic_map_alloc_div(bmap);
+		if (k < 0)
+			goto error;
+		isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
+		isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
+		if (isl_basic_map_add_div_constraints(bmap, k) < 0)
+			goto error;
+	}
+	k = isl_basic_map_alloc_equality(bmap);
+	if (k < 0)
+		goto error;
+	isl_int_neg(bmap->eq[k][pos], aff->el[0]);
+	isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
+	isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
+
+	isl_vec_free(aff);
+	isl_qpolynomial_free(qp);
+	bmap = isl_basic_map_finalize(bmap);
+	return bmap;
+error:
+	isl_vec_free(aff);
+	isl_qpolynomial_free(qp);
+	isl_basic_map_free(bmap);
+	return NULL;
+}