isl_qpolynomial_morph: properly handle denominators in morph
[platform/upstream/isl.git] / isl_polynomial.c
index 8acfa98..75908af 100644 (file)
@@ -9,16 +9,22 @@
  */
 
 #include <stdlib.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_constraint_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>
+#include <isl_aff_private.h>
+#include <isl_config.h>
 
 static unsigned pos(__isl_keep isl_dim *dim, enum isl_dim_type type)
 {
@@ -322,7 +328,7 @@ __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
        isl_assert(ctx, size >= 0, return NULL);
        rec = isl_calloc(ctx, struct isl_upoly_rec,
                        sizeof(struct isl_upoly_rec) +
-                       (size - 1) * sizeof(struct isl_upoly *));
+                       size * sizeof(struct isl_upoly *));
        if (!rec)
                return NULL;
 
@@ -470,8 +476,6 @@ error:
 
 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
 {
-       struct isl_upoly *dup;
-
        if (!up)
                return NULL;
 
@@ -674,7 +678,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;
 
@@ -687,12 +734,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;
@@ -701,7 +749,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);
@@ -709,7 +757,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;
        }
@@ -751,7 +799,7 @@ __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
 {
        struct isl_upoly_rec *rec1;
        struct isl_upoly_rec *rec2;
-       struct isl_upoly_rec *res;
+       struct isl_upoly_rec *res = NULL;
        int i, j;
        int size;
 
@@ -872,6 +920,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)
 {
@@ -960,10 +1033,9 @@ 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;
        struct isl_upoly_rec *rec;
        struct isl_upoly_cst *cst;
 
@@ -1004,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) {
@@ -1045,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;
@@ -1142,6 +1201,8 @@ static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
                    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++;
@@ -1169,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)
 {
@@ -1273,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;
 
@@ -1335,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,
@@ -1344,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;
 
-       qp->upoly = isl_upoly_neg(qp->upoly);
+       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_mul_isl_int(qp->upoly, v);
        if (!qp->upoly)
                goto error;
 
@@ -1361,6 +1419,12 @@ error:
        return NULL;
 }
 
+__isl_give isl_qpolynomial *isl_qpolynomial_scale(
+       __isl_take isl_qpolynomial *qp, isl_int v)
+{
+       return isl_qpolynomial_mul_isl_int(qp, v);
+}
+
 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
        __isl_take isl_qpolynomial *qp2)
 {
@@ -1389,28 +1453,56 @@ 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)
 {
+       if (!dim)
+               return NULL;
        return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_one(__isl_take isl_dim *dim)
 {
+       if (!dim)
+               return NULL;
        return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
 {
+       if (!dim)
+               return NULL;
        return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty(__isl_take isl_dim *dim)
 {
+       if (!dim)
+               return NULL;
        return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
 }
 
 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
 {
+       if (!dim)
+               return NULL;
        return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
 }
 
@@ -1420,6 +1512,9 @@ __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
        struct isl_qpolynomial *qp;
        struct isl_upoly_cst *cst;
 
+       if (!dim)
+               return NULL;
+
        qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
        if (!qp)
                return NULL;
@@ -1554,13 +1649,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)
@@ -1572,12 +1666,22 @@ error:
        return NULL;
 }
 
-int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
+int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
        __isl_keep isl_qpolynomial *qp2)
 {
+       int equal;
+
        if (!qp1 || !qp2)
                return -1;
 
+       equal = isl_dim_equal(qp1->dim, qp2->dim);
+       if (equal < 0 || !equal)
+               return equal;
+
+       equal = isl_mat_is_equal(qp1->div, qp2->div);
+       if (equal < 0 || !equal)
+               return equal;
+
        return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
 }
 
@@ -1611,7 +1715,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;
@@ -1621,7 +1725,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,
@@ -1636,12 +1740,81 @@ __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 struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
+       unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
+{
+       int i;
+       struct isl_upoly_rec *rec;
+       struct isl_upoly *base, *res;
+
+       if (!up)
+               return NULL;
+
+       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_assert(up->ctx, rec->n >= 1, goto error);
+
+       if (up->var >= first + n)
+               base = isl_upoly_var_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_var_pow(ctx, i, 1);
+               t = isl_upoly_mul(c, t);
+               up = isl_upoly_sum(up, t);
+       }
+
+       return up;
+}
+
 /* Remove common factor of non-constant terms and denominator.
  */
 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
@@ -1663,6 +1836,233 @@ static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
                            ctx->normalize_gcd);
 }
 
+/* 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;
+       int total;
+       int *reordering;
+
+       if (!qp || !s)
+               goto error;
+
+       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 or 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;
+       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)
 {
@@ -1684,10 +2084,8 @@ __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
        if (!qp)
                goto error;
 
-       for (i = 0; i < div->bmap->n_div; ++i) {
+       for (i = 0; i < div->bmap->n_div; ++i)
                isl_seq_cpy(qp->div->row[i], div->bmap->div[i], qp->div->n_col);
-               normalize_div(qp, i);
-       }
 
        for (i = 0; i < 1 + power; ++i) {
                rec->p[i] = isl_upoly_zero(div->ctx);
@@ -1700,6 +2098,8 @@ __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
 
        isl_div_free(div);
 
+       qp = reduce_divs(qp);
+
        return qp;
 error:
        isl_qpolynomial_free(qp);
@@ -1718,6 +2118,9 @@ __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
        struct isl_qpolynomial *qp;
        struct isl_upoly_cst *cst;
 
+       if (!dim)
+               return NULL;
+
        qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
        if (!qp)
                return NULL;
@@ -1787,7 +2190,7 @@ int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
        isl_assert(qp->dim->ctx, type == isl_dim_param ||
                                 type == isl_dim_set, return -1);
 
-       active = isl_calloc_array(set->ctx, int, isl_dim_total(qp->dim));
+       active = isl_calloc_array(qp->dim->ctx, int, isl_dim_total(qp->dim));
        if (set_active(qp, active) < 0)
                goto error;
 
@@ -1807,248 +2210,180 @@ error:
        return -1;
 }
 
-__isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
-       unsigned first, unsigned n)
-{
-       int i;
-       struct isl_upoly_rec *rec;
-
-       if (!up)
-               return NULL;
-       if (n == 0 || up->var < 0 || up->var < first)
-               return up;
-       if (up->var < first + n) {
-               up = replace_by_constant_term(up);
-               return isl_upoly_drop(up, first, n);
-       }
-       up = isl_upoly_cow(up);
-       if (!up)
-               return NULL;
-       up->var -= n;
-       rec = isl_upoly_as_rec(up);
-       if (!rec)
-               goto error;
-
-       for (i = 0; i < rec->n; ++i) {
-               rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
-               if (!rec->p[i])
-                       goto error;
-       }
-
-       return up;
-error:
-       isl_upoly_free(up);
-       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)
+/* Remove divs that do not appear in the quasi-polynomial, nor in any
+ * of the divs that do appear in the quasi-polynomial.
+ */
+static __isl_give isl_qpolynomial *remove_redundant_divs(
+       __isl_take isl_qpolynomial *qp)
 {
-       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;
-}
+       int i, j;
+       int d;
+       int len;
+       int skip;
+       int *active = NULL;
+       int *reordering = NULL;
+       int redundant = 0;
+       int n_div;
+       isl_ctx *ctx;
 
-__isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
-       __isl_take isl_qpolynomial *qp,
-       enum isl_dim_type type, unsigned first, unsigned n)
-{
        if (!qp)
                return NULL;
-       if (n == 0 && !isl_dim_get_tuple_name(qp->dim, type))
+       if (qp->div->n_row == 0)
                return qp;
 
-       qp = isl_qpolynomial_cow(qp);
-       if (!qp)
-               return NULL;
-
-       isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
-                       goto error);
-       isl_assert(qp->dim->ctx, type == isl_dim_param ||
-                                type == isl_dim_set, goto error);
-
-       qp->dim = isl_dim_drop(qp->dim, type, first, n);
-       if (!qp->dim)
+       d = isl_dim_total(qp->dim);
+       len = qp->div->n_col - 2;
+       ctx = isl_qpolynomial_get_ctx(qp);
+       active = isl_calloc_array(ctx, int, len);
+       if (!active)
                goto error;
 
-       if (type == isl_dim_set)
-               first += isl_dim_size(qp->dim, isl_dim_param);
-
-       qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
-       if (!qp->div)
+       if (up_set_active(qp->upoly, active, len) < 0)
                goto error;
 
-       qp->upoly = isl_upoly_drop(qp->upoly, first, n);
-       if (!qp->upoly)
-               goto error;
+       for (i = qp->div->n_row - 1; i >= 0; --i) {
+               if (!active[d + i]) {
+                       redundant = 1;
+                       continue;
+               }
+               for (j = 0; j < i; ++j) {
+                       if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
+                               continue;
+                       active[d + j] = 1;
+                       break;
+               }
+       }
 
-       return qp;
-error:
-       isl_qpolynomial_free(qp);
-       return NULL;
-}
+       if (!redundant) {
+               free(active);
+               return qp;
+       }
 
-__isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
-       unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
-{
-       int i;
-       struct isl_upoly_rec *rec;
-       struct isl_upoly *base, *res;
+       reordering = isl_alloc_array(qp->div->ctx, int, len);
+       if (!reordering)
+               goto error;
 
-       if (!up)
-               return NULL;
+       for (i = 0; i < d; ++i)
+               reordering[i] = i;
 
-       if (isl_upoly_is_cst(up))
-               return up;
+       skip = 0;
+       n_div = qp->div->n_row;
+       for (i = 0; i < n_div; ++i) {
+               if (!active[d + i]) {
+                       qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
+                       qp->div = isl_mat_drop_cols(qp->div,
+                                                   2 + d + i - skip, 1);
+                       skip++;
+               }
+               reordering[d + i] = d + i - skip;
+       }
 
-       if (up->var < first)
-               return up;
+       qp->upoly = reorder(qp->upoly, reordering);
 
-       rec = isl_upoly_as_rec(up);
-       if (!rec)
+       if (!qp->upoly || !qp->div)
                goto error;
 
-       isl_assert(up->ctx, rec->n >= 1, goto error);
-
-       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;
+       free(active);
+       free(reordering);
+
+       return qp;
 error:
-       isl_upoly_free(up);
+       free(active);
+       free(reordering);
+       isl_qpolynomial_free(qp);
        return NULL;
-}      
+}
 
-__isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
-       isl_int denom, unsigned len)
+__isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
+       unsigned first, unsigned n)
 {
        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;
+       struct isl_upoly_rec *rec;
 
-               if (isl_int_is_zero(f[1 + i]))
-                       continue;
+       if (!up)
+               return NULL;
+       if (n == 0 || up->var < 0 || up->var < first)
+               return up;
+       if (up->var < first + n) {
+               up = replace_by_constant_term(up);
+               return isl_upoly_drop(up, first, n);
+       }
+       up = isl_upoly_cow(up);
+       if (!up)
+               return NULL;
+       up->var -= n;
+       rec = isl_upoly_as_rec(up);
+       if (!rec)
+               goto error;
 
-               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);
+       for (i = 0; i < rec->n; ++i) {
+               rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
+               if (!rec->p[i])
+                       goto error;
        }
 
        return up;
+error:
+       isl_upoly_free(up);
+       return NULL;
 }
 
-/* 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_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
        __isl_take isl_qpolynomial *qp,
-       int div, __isl_take struct isl_upoly *s)
+       enum isl_dim_type type, unsigned pos, const char *s)
 {
-       int i;
-       int total;
-       int *reordering;
-
-       if (!qp || !s)
-               goto error;
-
        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;
-
-       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)
+__isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
+       __isl_take isl_qpolynomial *qp,
+       enum isl_dim_type type, unsigned first, unsigned n)
 {
-       int i, j;
-       int total;
-       struct isl_upoly *s;
+       if (!qp)
+               return NULL;
+       if (n == 0 && !isl_dim_is_named_or_nested(qp->dim, type))
+               return qp;
 
+       qp = isl_qpolynomial_cow(qp);
        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;
-       }
+       isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
+                       goto error);
+       isl_assert(qp->dim->ctx, type == isl_dim_param ||
+                                type == isl_dim_set, goto error);
+
+       qp->dim = isl_dim_drop(qp->dim, type, first, n);
+       if (!qp->dim)
+               goto error;
+
+       if (type == isl_dim_set)
+               first += isl_dim_size(qp->dim, isl_dim_param);
+
+       qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
+       if (!qp->div)
+               goto error;
+
+       qp->upoly = isl_upoly_drop(qp->upoly, first, n);
+       if (!qp->upoly)
+               goto error;
 
        return qp;
+error:
+       isl_qpolynomial_free(qp);
+       return NULL;
 }
 
-__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
+static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
        __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
 {
        int i, j, k;
@@ -2114,6 +2449,22 @@ error:
        return NULL;
 }
 
+/* Exploit the equalities in "eq" to simplify the quasi-polynomial.
+ */
+__isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
+       __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
+{
+       if (!qp || !eq)
+               goto error;
+       if (qp->div->n_row > 0)
+               eq = isl_basic_set_add(eq, isl_dim_set, qp->div->n_row);
+       return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
+error:
+       isl_basic_set_free(eq);
+       isl_qpolynomial_free(qp);
+       return NULL;
+}
+
 static __isl_give isl_basic_set *add_div_constraints(
        __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
 {
@@ -2163,7 +2514,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_gist(
        }
 
        aff = isl_set_affine_hull(context);
-       return isl_qpolynomial_substitute_equalities(qp, aff);
+       return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
 error:
        isl_qpolynomial_free(qp);
        isl_set_free(context);
@@ -2174,6 +2525,10 @@ error:
 #define PW isl_pw_qpolynomial
 #undef EL
 #define EL isl_qpolynomial
+#undef EL_IS_ZERO
+#define EL_IS_ZERO is_zero
+#undef ZERO
+#define ZERO zero
 #undef IS_ZERO
 #define IS_ZERO is_zero
 #undef FIELD
@@ -2198,7 +2553,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);
@@ -2210,7 +2565,6 @@ __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
 {
        int i, j, n;
        struct isl_pw_qpolynomial *res;
-       isl_set *set;
 
        if (!pwqp1 || !pwqp2)
                goto error;
@@ -2247,7 +2601,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;
                        }
@@ -2270,40 +2624,6 @@ error:
        return NULL;
 }
 
-__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
-       __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 = isl_qpolynomial_neg(pwqp->p[i].qp);
-               if (!pwqp->p[i].qp)
-                       goto error;
-       }
-
-       return pwqp;
-error:
-       isl_pw_qpolynomial_free(pwqp);
-       return NULL;
-}
-
-__isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
-       __isl_take isl_pw_qpolynomial *pwqp1,
-       __isl_take isl_pw_qpolynomial *pwqp2)
-{
-       return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
-}
-
 __isl_give struct isl_upoly *isl_upoly_eval(
        __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
 {
@@ -2483,7 +2803,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
        unsigned g_pos;
        int *exp;
 
-       if (n == 0)
+       if (n == 0 && !isl_dim_is_named_or_nested(qp->dim, type))
                return qp;
 
        qp = isl_qpolynomial_cow(qp);
@@ -2495,7 +2815,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
 
        g_pos = pos(qp->dim, type) + first;
 
-       qp->div = isl_mat_insert_cols(qp->div, 2 + g_pos, n);
+       qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
        if (!qp->div)
                goto error;
 
@@ -2638,39 +2958,47 @@ __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_dim *dim,
        return isl_qpolynomial_alloc(dim, 0, up);
 }
 
-__isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
-       __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
+__isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
 {
-       isl_int denom;
-       isl_dim *dim;
+       isl_ctx *ctx;
        struct isl_upoly *up;
        isl_qpolynomial *qp;
-       int sgn;
 
-       if (!c)
+       if (!aff)
                return NULL;
 
-       isl_int_init(denom);
-
-       isl_constraint_get_coefficient(c, type, pos, &denom);
-       isl_constraint_set_coefficient(c, type, pos, c->ctx->zero);
-       sgn = isl_int_sgn(denom);
-       isl_int_abs(denom, denom);
-       up = isl_upoly_from_affine(c->ctx, c->line[0], denom,
-                                       1 + isl_constraint_dim(c, isl_dim_all));
-       if (sgn < 0)
-               isl_int_neg(denom, denom);
-       isl_constraint_set_coefficient(c, type, pos, denom);
+       ctx = isl_aff_get_ctx(aff);
+       up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
+                                   aff->v->size - 1);
 
-       dim = isl_dim_copy(c->bmap->dim);
+       qp = isl_qpolynomial_alloc(isl_aff_get_dim(aff),
+                                   aff->ls->div->n_row, up);
+       if (!qp)
+               goto error;
 
-       isl_int_clear(denom);
-       isl_constraint_free(c);
+       isl_mat_free(qp->div);
+       qp->div = isl_mat_copy(aff->ls->div);
+       qp->div = isl_mat_cow(qp->div);
+       if (!qp->div)
+               goto error;
 
-       qp = isl_qpolynomial_alloc(dim, 0, up);
-       if (sgn > 0)
-               qp = isl_qpolynomial_neg(qp);
+       isl_aff_free(aff);
+       qp = reduce_divs(qp);
+       qp = remove_redundant_divs(qp);
        return qp;
+error:
+       isl_aff_free(aff);
+       return NULL;
+}
+
+__isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
+       __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
+{
+       isl_aff *aff;
+
+       aff = isl_constraint_get_bound(c, type, pos);
+       isl_constraint_free(c);
+       return isl_qpolynomial_from_aff(aff);
 }
 
 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
@@ -2908,7 +3236,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);
@@ -3226,7 +3554,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);
@@ -3412,10 +3740,8 @@ __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp
        int i;
        int n_sub;
        isl_ctx *ctx;
-       struct isl_upoly *up;
-       unsigned n_div;
        struct isl_upoly **subs;
-       isl_mat *mat;
+       isl_mat *mat, *diag;
 
        qp = isl_qpolynomial_cow(qp);
        if (!qp || !morph)
@@ -3437,7 +3763,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp
        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_pow(ctx, morph->inv->n_col - 1 + i, 1);
+                           isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
 
        qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
 
@@ -3445,8 +3771,10 @@ __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp
                isl_upoly_free(subs[i]);
        free(subs);
 
-       mat = isl_mat_diagonal(isl_mat_identity(ctx, 1), isl_mat_copy(morph->inv));
-       mat = isl_mat_diagonal(mat, isl_mat_identity(ctx, qp->div->n_row));
+       diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
+       mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
+       diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
+       mat = isl_mat_diagonal(mat, diag);
        qp->div = isl_mat_product(qp->div, mat);
        isl_dim_free(qp->dim);
        qp->dim = isl_dim_copy(morph->ran->dim);
@@ -3603,6 +3931,33 @@ error:
        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;
@@ -3935,7 +4290,7 @@ __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
        if (!bset)
                return NULL;
 
-       if (isl_basic_set_fast_is_empty(bset))
+       if (isl_basic_set_plain_is_empty(bset))
                return constant_on_domain(bset, 0);
 
        orig_nvar = isl_basic_set_dim(bset, isl_dim_set);
@@ -4236,3 +4591,55 @@ 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;
+}