Merge branch 'maint'
[platform/upstream/isl.git] / isl_polynomial.c
index f5f1828..7eae4f8 100644 (file)
@@ -9,6 +9,8 @@
  */
 
 #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_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 +327,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 +475,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 +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;
 
@@ -687,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;
@@ -701,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);
@@ -709,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;
        }
@@ -751,7 +798,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;
 
@@ -988,7 +1035,6 @@ void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
 __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;
 
@@ -1070,19 +1116,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;
@@ -1196,50 +1229,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)
 {
@@ -1300,7 +1289,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;
 
@@ -1362,7 +1351,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,
@@ -1371,14 +1362,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;
 
@@ -1436,26 +1466,36 @@ error:
 
 __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));
 }
 
@@ -1465,6 +1505,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;
@@ -1599,13 +1642,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)
@@ -1687,6 +1729,75 @@ 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)
+{
+       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)
@@ -1708,6 +1819,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)
 {
@@ -1729,10 +2067,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);
@@ -1745,7 +2081,7 @@ __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
 
        isl_div_free(div);
 
-       qp = sort_divs(qp);
+       qp = reduce_divs(qp);
 
        return qp;
 error:
@@ -1834,24 +2170,108 @@ 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));
-       if (set_active(qp, active) < 0)
+       active = isl_calloc_array(qp->dim->ctx, int, isl_dim_total(qp->dim));
+       if (set_active(qp, active) < 0)
+               goto error;
+
+       if (type == isl_dim_set)
+               first += isl_dim_size(qp->dim, isl_dim_param);
+       for (i = 0; i < n; ++i)
+               if (active[first + i]) {
+                       involves = 1;
+                       break;
+               }
+
+       free(active);
+
+       return involves;
+error:
+       free(active);
+       return -1;
+}
+
+/* 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)
+{
+       int i, j;
+       int d;
+       int len;
+       int skip;
+       int *active = NULL;
+       int *reordering = NULL;
+       int redundant = 0;
+       int n_div;
+       isl_ctx *ctx;
+
+       if (!qp)
+               return NULL;
+       if (qp->div->n_row == 0)
+               return qp;
+
+       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 (up_set_active(qp->upoly, active, len) < 0)
                goto error;
 
-       if (type == isl_dim_set)
-               first += isl_dim_size(qp->dim, isl_dim_param);
-       for (i = 0; i < n; ++i)
-               if (active[first + i]) {
-                       involves = 1;
+       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;
                }
+       }
+
+       if (!redundant) {
+               free(active);
+               return qp;
+       }
+
+       reordering = isl_alloc_array(qp->div->ctx, int, len);
+       if (!reordering)
+               goto error;
+
+       for (i = 0; i < d; ++i)
+               reordering[i] = i;
+
+       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;
+       }
+
+       qp->upoly = reorder(qp->upoly, reordering);
+
+       if (!qp->upoly || !qp->div)
+               goto error;
 
        free(active);
+       free(reordering);
 
-       return involves;
+       return qp;
 error:
        free(active);
-       return -1;
+       free(reordering);
+       isl_qpolynomial_free(qp);
+       return NULL;
 }
 
 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
@@ -1943,158 +2363,6 @@ 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)
-{
-       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;
-}
-
-/* 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;
-}
-
 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
        __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
 {
@@ -2245,7 +2513,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);
@@ -2257,7 +2525,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;
@@ -2294,7 +2561,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;
                        }
@@ -2685,6 +2952,39 @@ __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_aff(__isl_take isl_aff *aff)
+{
+       isl_ctx *ctx;
+       struct isl_upoly *up;
+       isl_qpolynomial *qp;
+
+       if (!aff)
+               return NULL;
+
+       ctx = isl_aff_get_ctx(aff);
+       up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
+                                   aff->v->size - 1);
+
+       qp = isl_qpolynomial_alloc(isl_aff_get_dim(aff),
+                                   aff->ls->div->n_row, up);
+       if (!qp)
+               goto error;
+
+       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;
+
+       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)
 {
@@ -3459,8 +3759,6 @@ __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;
 
@@ -3650,6 +3948,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;
@@ -3982,7 +4307,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);
@@ -4283,3 +4608,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;
+}