isl_basic_{map,set}: explicitly store constraints defining div

These constraints are redundant, but we need then while performing
various operations, so it is easier to keep them available explicitly.

diff --git a/include/isl_int.h b/include/isl_int.h
index ebc28c5..731e240 100644 (file)
--- a/include/isl_int.h
+++ b/include/isl_int.h
@@ -23,6 +23,8 @@ typedef mpz_t isl_int;
 #define isl_int_abs(r,i)       mpz_abs(r,i)
 #define isl_int_neg(r,i)       mpz_neg(r,i)
 #define isl_int_swap(i,j)      mpz_swap(i,j)
+#define isl_int_swap_or_set(i,j)       mpz_swap(i,j)
+#define isl_int_add_ui(r,i,j)  mpz_add_ui(r,i,j)
 #define isl_int_sub_ui(r,i,j)  mpz_sub_ui(r,i,j)
 
 #define isl_int_add(r,i,j)     mpz_add(r,i,j)
@@ -51,6 +53,8 @@ typedef mpz_t isl_int;
 #define isl_int_eq(i,j)                (mpz_cmp(i,j) == 0)
 #define isl_int_ne(i,j)                (mpz_cmp(i,j) != 0)
 #define isl_int_lt(i,j)                (mpz_cmp(i,j) < 0)
+#define isl_int_abs_eq(i,j)    (mpz_cmpabs(i,j) == 0)
+#define isl_int_abs_ne(i,j)    (mpz_cmpabs(i,j) != 0)
 #define isl_int_abs_lt(i,j)    (mpz_cmpabs(i,j) < 0)
 
 

diff --git a/include/isl_seq.h b/include/isl_seq.h
index 98864b7..120182b 100644 (file)
--- a/include/isl_seq.h
+++ b/include/isl_seq.h
@@ -17,6 +17,7 @@ void isl_seq_gcd(isl_int *p, unsigned len, isl_int *gcd);
 int isl_seq_first_non_zero(isl_int *p, unsigned len);
 int isl_seq_abs_min_non_zero(isl_int *p, unsigned len);
 int isl_seq_eq(isl_int *p1, isl_int *p2, unsigned len);
+int isl_seq_is_neg(isl_int *p1, isl_int *p2, unsigned len);
 
 u_int32_t isl_seq_hash(isl_int *p, unsigned len, unsigned bits);
 

diff --git a/isl_map.c b/isl_map.c
index 2df3161..66cba37 100644 (file)
--- a/isl_map.c
+++ b/isl_map.c
@@ -590,8 +590,14 @@ static struct isl_basic_map *isl_basic_map_drop_div(struct isl_ctx *ctx,
        for (i = 0; i < bmap->n_eq; ++i)
                constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
 
-       for (i = 0; i < bmap->n_ineq; ++i)
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               if (!isl_int_is_zero(bmap->ineq[i][pos])) {
+                       isl_basic_map_drop_inequality(ctx, bmap, i);
+                       --i;
+                       continue;
+               }
                constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
+       }
 
        for (i = 0; i < bmap->n_div; ++i)
                constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
@@ -2209,6 +2215,14 @@ error:
        return NULL;
 }
 
+/* If the only constraints a div d=floor(f/m)
+ * appears in are its two defining constraints
+ *
+ *     f - m d >=0
+ *     -(f - (m - 1)) + m d >= 0
+ *
+ * then it can safely be removed.
+ */
 static int div_is_redundant(struct isl_ctx *ctx, struct isl_basic_map *bmap,
                                int div)
 {
@@ -2219,9 +2233,32 @@ static int div_is_redundant(struct isl_ctx *ctx, struct isl_basic_map *bmap,
                if (!isl_int_is_zero(bmap->eq[i][pos]))
                        return 0;
 
-       for (i = 0; i < bmap->n_ineq; ++i)
-               if (!isl_int_is_zero(bmap->ineq[i][pos]))
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               if (isl_int_is_zero(bmap->ineq[i][pos]))
+                       continue;
+               if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
+                       int neg;
+                       isl_int_sub(bmap->div[div][1],
+                                       bmap->div[div][1], bmap->div[div][0]);
+                       isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
+                       neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
+                       isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
+                       isl_int_add(bmap->div[div][1],
+                                       bmap->div[div][1], bmap->div[div][0]);
+                       if (!neg)
+                               return 0;
+                       if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
+                                                   bmap->n_div-div-1) != -1)
+                               return 0;
+               } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
+                       if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
+                               return 0;
+                       if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
+                                                   bmap->n_div-div-1) != -1)
+                               return 0;
+               } else
                        return 0;
+       }
 
        for (i = 0; i < bmap->n_div; ++i)
                if (!isl_int_is_zero(bmap->div[i][1+pos]))
@@ -2549,6 +2586,37 @@ static int find_div(struct isl_basic_map *dst,
        return -1;
 }
 
+/* For a div d = floor(f/m), add the constraints
+ *
+ *             f - m d >= 0
+ *             -(f-(n-1)) + m d >= 0
+ *
+ * Note that the second constraint is the negation of
+ *
+ *             f - m d >= n
+ */
+static int add_div_constraints(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned div)
+{
+       int i, j;
+       unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       i = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (i < 0)
+               return -1;
+       isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
+       isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
+
+       j = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (j < 0)
+               return -1;
+       isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
+       isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
+       isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
+       return j;
+}
+
 struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
                struct isl_basic_map *dst, struct isl_basic_map *src)
 {
@@ -2560,7 +2628,7 @@ struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
 
        src = order_divs(ctx, src);
        dst = isl_basic_map_extend(ctx, dst, dst->nparam, dst->n_in,
-                       dst->n_out, src->n_div, 0, 0);
+                       dst->n_out, src->n_div, 0, 2 * src->n_div);
        if (!dst)
                return NULL;
        for (i = 0; i < src->n_div; ++i) {
@@ -2570,6 +2638,8 @@ struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
                        if (j < 0)
                                goto error;
                        isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
+                       if (add_div_constraints(ctx, dst, j) < 0)
+                               goto error;
                }
                if (j != i)
                        swap_div(dst, i, j);

diff --git a/isl_map_piplib.c b/isl_map_piplib.c
index f5741c0..83617e9 100644 (file)
--- a/isl_map_piplib.c
+++ b/isl_map_piplib.c
@@ -45,6 +45,37 @@ static int add_inequality(struct isl_ctx *ctx,
        return i;
 }
 
+/* For a div d = floor(f/m), add the constraints
+ *
+ *             f - m d >= 0
+ *             -(f-(n-1)) + m d >= 0
+ *
+ * Note that the second constraint is the negation of
+ *
+ *             f - m d >= n
+ */
+static int add_div_constraints(struct isl_ctx *ctx,
+       struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
+{
+       int i, j;
+       unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+
+       i = add_inequality(ctx, bmap, pos, p->vector);
+       if (i < 0)
+               return -1;
+       copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
+       isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
+
+       j = isl_basic_map_alloc_inequality(ctx, bmap);
+       if (j < 0)
+               return -1;
+       isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
+       isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
+       isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
+       return j;
+}
+
 static int add_equality(struct isl_ctx *ctx,
                   struct isl_basic_map *bmap, int *pos,
                   unsigned var, PipVector *vec)
@@ -69,7 +100,8 @@ static int find_div(struct isl_ctx *ctx,
        int i, j;
 
        i = isl_basic_map_alloc_div(ctx, bmap);
-       assert(i != -1);
+       if (i < 0)
+               return -1;
 
        copy_constraint_from(bmap->div[i]+1, p->vector,
            bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, pos);
@@ -82,6 +114,9 @@ static int find_div(struct isl_ctx *ctx,
                        return j;
                }
 
+       if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
+               return -1;
+
        return i;
 }
 
@@ -141,7 +176,6 @@ static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
 
        for (p = q->newparm; p; p = p->next) {
                int pos;
-               PipVector *vec = p->vector;
                unsigned pip_param = bmap->nparam + bmap->n_in;
 
                pos = find_div(data->ctx, bmap, data->pos, p);
@@ -191,6 +225,9 @@ static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
                *data->rest = isl_set_add(data->ctx, *data->rest, bset);
        }
 
+       if (isl_basic_map_free_inequality(data->ctx, bmap,
+                                          2*(bmap->n_div - old_n_div)))
+               goto error;
        if (isl_basic_map_free_div(data->ctx, bmap,
                                           bmap->n_div - old_n_div))
                goto error;
@@ -249,7 +286,7 @@ static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
        data.bmap = isl_basic_map_from_basic_set(ctx, context,
                                context->dim, 0);
        data.bmap = isl_basic_map_extend(data.ctx, data.bmap,
-                       nparam, data.bmap->n_in, keep, nexist, keep, max_depth);
+           nparam, data.bmap->n_in, keep, nexist, keep, max_depth+2*nexist);
        if (!data.bmap)
                goto error;
 
@@ -316,7 +353,7 @@ PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
                                + bmap->n_div - pip_param;
        unsigned pip_dim = pip_var - bmap->n_div;
 
-       nrow = extra_front + bmap->n_eq + bmap->n_ineq + 2*bmap->n_div;
+       nrow = extra_front + bmap->n_eq + bmap->n_ineq;
        ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
        M = pip_matrix_alloc(nrow, ncol);
        if (!M)
@@ -334,29 +371,6 @@ PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
                copy_constraint_to(M->p[off+i], bmap->ineq[i],
                                   pip_param, pip_var, extra_front, extra_back);
        }
-       off += bmap->n_ineq;
-       for (i = 0; i < bmap->n_div; ++i) {
-               unsigned total = bmap->n_in + bmap->n_out
-                                   + bmap->n_div + bmap->nparam + extra_back;
-               if (isl_int_is_zero(bmap->div[i][0]))
-                       continue;
-               entier_set_si(M->p[off+2*i][0], 1);
-               copy_constraint_to(M->p[off+2*i], bmap->div[i]+1,
-                                  pip_param, pip_var, extra_front, extra_back);
-               copy_values_to(M->p[off+2*i]+1+extra_front+pip_dim+i,
-                               bmap->div[i], 1);
-               entier_oppose(M->p[off+2*i][1+extra_front+pip_dim+i],
-                             M->p[off+2*i][1+extra_front+pip_dim+i]);
-
-               entier_set_si(M->p[off+2*i+1][0], 1);
-               Vector_Oppose(M->p[off+2*i]+1+extra_front,
-                             M->p[off+2*i+1]+1+extra_front, total+1);
-               entier_addto(M->p[off+2*i+1][1+extra_front+total],
-                            M->p[off+2*i+1][1+extra_front+total],
-                            M->p[off+2*i+1][1+extra_front+pip_dim+i]);
-               entier_decrement(M->p[off+2*i+1][1+extra_front+total],
-                                M->p[off+2*i+1][1+extra_front+total]);
-       }
        return M;
 }
 

diff --git a/isl_map_polylib.c b/isl_map_polylib.c
index 58848ac..3502f82 100644 (file)
--- a/isl_map_polylib.c
+++ b/isl_map_polylib.c
@@ -163,7 +163,7 @@ Polyhedron *isl_basic_map_to_polylib(struct isl_ctx *ctx,
        Polyhedron *P;
        unsigned off;
 
-       M = Matrix_Alloc(bmap->n_eq + bmap->n_ineq + 2*bmap->n_div,
+       M = Matrix_Alloc(bmap->n_eq + bmap->n_ineq,
                 1 + bmap->n_in + bmap->n_out + bmap->n_div + bmap->nparam + 1);
        for (i = 0; i < bmap->n_eq; ++i) {
                value_set_si(M->p[i][0], 0);
@@ -176,26 +176,6 @@ Polyhedron *isl_basic_map_to_polylib(struct isl_ctx *ctx,
                copy_constraint_to(M->p[off+i], bmap->ineq[i],
                           bmap->nparam, bmap->n_in + bmap->n_out, bmap->n_div);
        }
-       off += bmap->n_ineq;
-       for (i = 0; i < bmap->n_div; ++i) {
-               unsigned total = bmap->n_in+bmap->n_out+bmap->n_div+bmap->nparam;
-               if (isl_int_is_zero(bmap->div[i][0]))
-                       continue;
-               value_set_si(M->p[off+2*i][0], 1);
-               copy_constraint_to(M->p[off+2*i], bmap->div[i]+1,
-                          bmap->nparam, bmap->n_in + bmap->n_out, bmap->n_div);
-               copy_values_to(M->p[off+2*i]+1+bmap->n_in+bmap->n_out+i,
-                               bmap->div[i], 1);
-               value_oppose(M->p[off+2*i][1+bmap->n_in+bmap->n_out+i],
-                            M->p[off+2*i][1+bmap->n_in+bmap->n_out+i]);
-
-               value_set_si(M->p[off+2*i+1][0], 1);
-               Vector_Oppose(M->p[off+2*i]+1, M->p[off+2*i+1]+1, total+1);
-               value_addto(M->p[off+2*i+1][1+total], M->p[off+2*i+1][1+total],
-                           M->p[off+2*i+1][1+bmap->n_in+bmap->n_out+i]);
-               value_decrement(M->p[off+2*i+1][1+total],
-                               M->p[off+2*i+1][1+total]);
-       }
        P = Constraints2Polyhedron(M, ctx->MaxRays);
        Matrix_Free(M);
 

diff --git a/isl_seq.c b/isl_seq.c
index 621fe85..5d02804 100644 (file)
--- a/isl_seq.c
+++ b/isl_seq.c
@@ -90,6 +90,21 @@ int isl_seq_eq(isl_int *p1, isl_int *p2, unsigned len)
        return 1;
 }
 
+int isl_seq_is_neg(isl_int *p1, isl_int *p2, unsigned len)
+{
+       int i;
+
+       for (i = 0; i < len; ++i) {
+               if (isl_int_abs_ne(p1[i], p2[i]))
+                       return 0;
+               if (isl_int_is_zero(p1[i]))
+                       continue;
+               if (isl_int_eq(p1[i], p2[i]))
+                       return 0;
+       }
+       return 1;
+}
+
 int isl_seq_first_non_zero(isl_int *p, unsigned len)
 {
        int i;