extract isl_map_align_divs from isl_affine_hull.c
[platform/upstream/isl.git] / isl_map.c
index 0a8fe3b..033a1d3 100644 (file)
--- a/isl_map.c
+++ b/isl_map.c
@@ -7,6 +7,8 @@
 #include "isl_map.h"
 #include "isl_map_private.h"
 #include "isl_map_piplib.h"
+#include "isl_sample.h"
+#include "isl_vec.h"
 
 static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
                struct isl_basic_map *bmap,
@@ -17,7 +19,7 @@ static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
        size_t row_size = 1 + nparam + n_in + n_out + extra;
 
        bmap->block = isl_blk_alloc(ctx, (n_eq + n_ineq) * row_size);
-       if (!bmap->block.data) {
+       if (isl_blk_is_error(bmap->block)) {
                free(bmap);
                return NULL;
        }
@@ -30,12 +32,11 @@ static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
        }
 
        if (extra == 0) {
-               bmap->block2.size = 0;
-               bmap->block2.data = NULL;
+               bmap->block2 = isl_blk_empty();
                bmap->div = NULL;
        } else {
                bmap->block2 = isl_blk_alloc(ctx, extra * (1 + row_size));
-               if (!bmap->block2.data) {
+               if (isl_blk_is_error(bmap->block2)) {
                        free(bmap->eq);
                        isl_blk_free(ctx, bmap->block);
                        free(bmap);
@@ -69,6 +70,7 @@ static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
        bmap->n_eq = 0;
        bmap->n_ineq = 0;
        bmap->n_div = 0;
+       bmap->sample = NULL;
 
        return bmap;
 }
@@ -130,6 +132,7 @@ struct isl_basic_map *isl_basic_map_dup(struct isl_ctx *ctx,
        if (!dup)
                return NULL;
        dup_constraints(ctx, dup, bmap);
+       dup->sample = isl_vec_copy(ctx, bmap->sample);
        return dup;
 }
 
@@ -198,6 +201,7 @@ void isl_basic_map_free(struct isl_ctx *ctx, struct isl_basic_map *bmap)
        isl_blk_free(ctx, bmap->block2);
        free(bmap->eq);
        isl_blk_free(ctx, bmap->block);
+       isl_vec_free(ctx, bmap->sample);
        free(bmap);
 }
 
@@ -227,6 +231,12 @@ int isl_basic_map_alloc_equality(struct isl_ctx *ctx,
        return bmap->n_eq++;
 }
 
+int isl_basic_set_alloc_equality(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       return isl_basic_map_alloc_equality(ctx, (struct isl_basic_map *)bset);
+}
+
 int isl_basic_map_free_equality(struct isl_ctx *ctx,
                struct isl_basic_map *bmap, unsigned n)
 {
@@ -250,14 +260,35 @@ int isl_basic_map_drop_equality(struct isl_ctx *ctx,
        return 0;
 }
 
+void isl_basic_map_inequality_to_equality(struct isl_ctx *ctx,
+               struct isl_basic_map *bmap, unsigned pos)
+{
+       isl_int *t;
+
+       t = bmap->ineq[pos];
+       bmap->ineq[pos] = bmap->ineq[0];
+       bmap->ineq[0] = bmap->eq[bmap->n_eq];
+       bmap->eq[bmap->n_eq] = t;
+       bmap->n_eq++;
+       bmap->n_ineq--;
+       bmap->ineq++;
+}
+
 int isl_basic_map_alloc_inequality(struct isl_ctx *ctx,
                struct isl_basic_map *bmap)
 {
        isl_assert(ctx, (bmap->ineq - bmap->eq) + bmap->n_ineq < bmap->c_size,
                        return -1);
+       F_CLR(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
        return bmap->n_ineq++;
 }
 
+int isl_basic_set_alloc_inequality(struct isl_ctx *ctx,
+               struct isl_basic_set *bset)
+{
+       return isl_basic_map_alloc_inequality(ctx, (struct isl_basic_map *)bset);
+}
+
 int isl_basic_map_free_inequality(struct isl_ctx *ctx,
                struct isl_basic_map *bmap, unsigned n)
 {
@@ -307,15 +338,31 @@ static void copy_constraint(struct isl_basic_map *dst_map, isl_int *dst,
 {
        isl_int_set(dst[0], src[0]);
        isl_seq_cpy(dst+1, src+1, isl_min(dst_map->nparam, src_map->nparam));
+       if (dst_map->nparam > src_map->nparam)
+               isl_seq_clr(dst+1+src_map->nparam,
+                               dst_map->nparam - src_map->nparam);
+       isl_seq_clr(dst+1+dst_map->nparam, in_off);
        isl_seq_cpy(dst+1+dst_map->nparam+in_off,
                    src+1+src_map->nparam,
                    isl_min(dst_map->n_in-in_off, src_map->n_in));
+       if (dst_map->n_in-in_off > src_map->n_in)
+               isl_seq_clr(dst+1+dst_map->nparam+in_off+src_map->n_in,
+                               dst_map->n_in - in_off - src_map->n_in);
+       isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in, out_off);
        isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+out_off,
                    src+1+src_map->nparam+src_map->n_in,
                    isl_min(dst_map->n_out-out_off, src_map->n_out));
+       if (dst_map->n_out-out_off > src_map->n_out)
+               isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+out_off+src_map->n_out,
+                               dst_map->n_out - out_off - src_map->n_out);
+       isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out, div_off);
        isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+div_off,
                    src+1+src_map->nparam+src_map->n_in+src_map->n_out,
                    isl_min(dst_map->extra-div_off, src_map->extra));
+       if (dst_map->extra-div_off > src_map->extra)
+               isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+
+                               div_off+src_map->extra,
+                               dst_map->extra - div_off - src_map->extra);
 }
 
 static void copy_div(struct isl_basic_map *dst_map, isl_int *dst,
@@ -421,7 +468,7 @@ struct isl_basic_set *isl_basic_set_extend(struct isl_ctx *ctx,
                                        nparam, 0, dim, extra, n_eq, n_ineq);
 }
 
-static struct isl_basic_set *isl_basic_set_cow(struct isl_ctx *ctx,
+struct isl_basic_set *isl_basic_set_cow(struct isl_ctx *ctx,
                struct isl_basic_set *bset)
 {
        return (struct isl_basic_set *)
@@ -444,6 +491,9 @@ struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
 
 static struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
 {
+       if (!set)
+               return NULL;
+
        if (set->ref == 1)
                return set;
        set->ref--;
@@ -452,6 +502,9 @@ static struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
 
 struct isl_map *isl_map_cow(struct isl_ctx *ctx, struct isl_map *map)
 {
+       if (!map)
+               return NULL;
+
        if (map->ref == 1)
                return map;
        map->ref--;
@@ -485,7 +538,7 @@ struct isl_basic_set *isl_basic_set_swap_vars(struct isl_ctx *ctx,
                return NULL;
 
        blk = isl_blk_alloc(ctx, bset->dim);
-       if (!blk.data)
+       if (isl_blk_is_error(blk))
                goto error;
 
        for (i = 0; i < bset->n_eq; ++i)
@@ -653,7 +706,7 @@ struct isl_basic_map *isl_basic_map_set_to_empty(
        unsigned total;
        if (!bmap)
                goto error;
-       total = bmap->nparam + bmap->n_in + bmap->n_out;
+       total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
        isl_basic_map_free_div(ctx, bmap, bmap->n_div);
        isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
        if (bmap->n_eq > 0)
@@ -672,6 +725,13 @@ error:
        return NULL;
 }
 
+struct isl_basic_set *isl_basic_set_set_to_empty(
+               struct isl_ctx *ctx, struct isl_basic_set *bset)
+{
+       return (struct isl_basic_set *)
+               isl_basic_map_set_to_empty(ctx, (struct isl_basic_map *)bset);
+}
+
 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
 {
        isl_int *t = bmap->eq[a];
@@ -786,6 +846,13 @@ struct isl_basic_map *isl_basic_map_gauss(struct isl_ctx *ctx,
        return bmap;
 }
 
+struct isl_basic_set *isl_basic_set_gauss(struct isl_ctx *ctx,
+       struct isl_basic_set *bset, int *progress)
+{
+       return (struct isl_basic_set*)isl_basic_map_gauss(ctx,
+                       (struct isl_basic_map *)bset, progress);
+}
+
 static unsigned int round_up(unsigned int v)
 {
        int old_v = v;
@@ -797,6 +864,20 @@ static unsigned int round_up(unsigned int v)
        return old_v << 1;
 }
 
+static int hash_index(int *index, unsigned int size, int bits,
+                       struct isl_basic_map *bmap, int k)
+{
+       int h;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
+       for (h = hash; index[h]; h = (h+1) % size)
+               if (k != index[h]-1 &&
+                   isl_seq_eq(bmap->ineq[k]+1,
+                              bmap->ineq[index[h]-1]+1, total))
+                       break;
+       return h;
+}
+
 static struct isl_basic_map *remove_duplicate_constraints(
        struct isl_ctx *ctx,
        struct isl_basic_map *bmap, int *progress)
@@ -806,11 +887,12 @@ static struct isl_basic_map *remove_duplicate_constraints(
        int k, l, h;
        int bits;
        unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       isl_int sum;
 
        if (bmap->n_ineq <= 1)
                return bmap;
 
-       size = round_up(4 * bmap->n_ineq / 3 - 1);
+       size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
        bits = ffs(size) - 1;
        index = isl_alloc_array(ctx, int, size);
        memset(index, 0, size * sizeof(int));
@@ -819,11 +901,7 @@ static struct isl_basic_map *remove_duplicate_constraints(
 
        index[isl_seq_hash(bmap->ineq[0]+1, total, bits)] = 1;
        for (k = 1; k < bmap->n_ineq; ++k) {
-               u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
-               for (h = hash; index[h]; h = (h+1) % size)
-                       if (isl_seq_eq(bmap->ineq[k]+1,
-                                      bmap->ineq[index[h]-1]+1, total))
-                               break;
+               h = hash_index(index, size, bits, bmap, k);
                if (!index[h]) {
                        index[h] = k+1;
                        continue;
@@ -835,6 +913,30 @@ static struct isl_basic_map *remove_duplicate_constraints(
                isl_basic_map_drop_inequality(ctx, bmap, k);
                --k;
        }
+       isl_int_init(sum);
+       for (k = 0; k < bmap->n_ineq-1; ++k) {
+               isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
+               h = hash_index(index, size, bits, bmap, k);
+               isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
+               if (!index[h])
+                       continue;
+               l = index[h] - 1;
+               isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
+               if (isl_int_is_pos(sum))
+                       continue;
+               if (isl_int_is_zero(sum)) {
+                       /* We need to break out of the loop after these
+                        * changes since the contents of the hash
+                        * will no longer be valid.
+                        * Plus, we probably we want to regauss first.
+                        */
+                       isl_basic_map_drop_inequality(ctx, bmap, l);
+                       isl_basic_map_inequality_to_equality(ctx, bmap, k);
+               } else
+                       bmap = isl_basic_map_set_to_empty(ctx, bmap);
+               break;
+       }
+       isl_int_clear(sum);
 
        free(index);
        return bmap;
@@ -867,7 +969,7 @@ static struct isl_basic_map *remove_duplicate_divs(struct isl_ctx *ctx,
        if (!index)
                return bmap;
        eq = isl_blk_alloc(ctx, 1+total);
-       if (!eq.data)
+       if (isl_blk_is_error(eq))
                goto out;
 
        isl_seq_clr(eq.data, 1+total);
@@ -1050,6 +1152,8 @@ static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
 static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
                                const char *op, FILE *out, int indent)
 {
+       int i;
+
        fprintf(out, "%*s", indent, "");
 
        dump_constraint_sign(bmap, c, 1, out);
@@ -1057,6 +1161,13 @@ static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
        dump_constraint_sign(bmap, c, -1, out);
 
        fprintf(out, "\n");
+
+       for (i = bmap->n_div; i < bmap->extra; ++i) {
+               if (isl_int_is_zero(c[1+bmap->nparam+bmap->n_in+bmap->n_out+i]))
+                       continue;
+               fprintf(out, "%*s", indent, "");
+               fprintf(out, "ERROR: unused div coefficient not zero\n");
+       }
 }
 
 static void dump_constraints(struct isl_basic_map *bmap,
@@ -1105,6 +1216,11 @@ static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
 void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
                                FILE *out, int indent)
 {
+       if (!bset) {
+               fprintf(out, "null basic set\n");
+               return;
+       }
+
        fprintf(out, "%*s", indent, "");
        fprintf(out, "nparam: %d, dim: %d, extra: %d\n",
                        bset->nparam, bset->dim, bset->extra);
@@ -1114,6 +1230,11 @@ void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
 void isl_basic_map_dump(struct isl_ctx *ctx, struct isl_basic_map *bmap,
                                FILE *out, int indent)
 {
+       if (!bmap) {
+               fprintf(out, "null basic map\n");
+               return;
+       }
+
        fprintf(out, "%*s", indent, "");
        fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d\n",
                bmap->ref,
@@ -1251,6 +1372,11 @@ void isl_map_dump(struct isl_ctx *ctx, struct isl_map *map, FILE *out,
 {
        int i;
 
+       if (!map) {
+               fprintf(out, "null map\n");
+               return;
+       }
+
        fprintf(out, "%*s", indent, "");
        fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d\n",
                        map->ref, map->n, map->nparam, map->n_in, map->n_out);
@@ -1282,6 +1408,7 @@ struct isl_basic_map *isl_basic_map_intersect_domain(
                                                bset->dim, 0);
        bmap = add_constraints(ctx, bmap, bmap_domain, 0, 0);
 
+       bmap = isl_basic_map_simplify(ctx, bmap);
        return isl_basic_map_finalize(ctx, bmap);
 error:
        isl_basic_map_free(ctx, bmap);
@@ -1310,6 +1437,7 @@ struct isl_basic_map *isl_basic_map_intersect_range(
                                                0, bset->dim);
        bmap = add_constraints(ctx, bmap, bmap_range, 0, 0);
 
+       bmap = isl_basic_map_simplify(ctx, bmap);
        return isl_basic_map_finalize(ctx, bmap);
 error:
        isl_basic_map_free(ctx, bmap);
@@ -1321,6 +1449,9 @@ struct isl_basic_map *isl_basic_map_intersect(
                struct isl_ctx *ctx, struct isl_basic_map *bmap1,
                struct isl_basic_map *bmap2)
 {
+       if (!bmap1 || !bmap2)
+               goto error;
+
        isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
        isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
        isl_assert(ctx, bmap1->n_out == bmap2->n_out, goto error);
@@ -1332,6 +1463,7 @@ struct isl_basic_map *isl_basic_map_intersect(
                goto error;
        bmap1 = add_constraints(ctx, bmap1, bmap2, 0, 0);
 
+       bmap1 = isl_basic_map_simplify(ctx, bmap1);
        return isl_basic_map_finalize(ctx, bmap1);
 error:
        isl_basic_map_free(ctx, bmap1);
@@ -1480,6 +1612,7 @@ struct isl_basic_map *isl_basic_map_apply_range(
        if (!bmap1)
                goto error;
        bmap1 = add_constraints(ctx, bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
+       bmap1 = isl_basic_map_simplify(ctx, bmap1);
        bset = isl_basic_set_from_basic_map(ctx, bmap1);
        bset = isl_basic_set_project_out(ctx, bset,
                                                bset->dim - (n_in + n_out), 0);
@@ -1643,6 +1776,26 @@ error:
        return NULL;
 }
 
+static struct isl_basic_set *isl_basic_map_underlying_set(
+               struct isl_ctx *ctx, struct isl_basic_map *bmap)
+{
+       if (!bmap)
+               goto error;
+       if (bmap->nparam == 0 && bmap->n_in == 0 && bmap->n_div == 0)
+               return (struct isl_basic_set *)bmap;
+       bmap = isl_basic_map_cow(ctx, bmap);
+       if (!bmap)
+               goto error;
+       bmap->n_out += bmap->nparam + bmap->n_in + bmap->n_div;
+       bmap->nparam = 0;
+       bmap->n_in = 0;
+       bmap->n_div = 0;
+       bmap = isl_basic_map_finalize(ctx, bmap);
+       return (struct isl_basic_set *)bmap;
+error:
+       return NULL;
+}
+
 struct isl_basic_set *isl_basic_map_domain(struct isl_ctx *ctx,
                struct isl_basic_map *bmap)
 {
@@ -1849,8 +2002,11 @@ struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_ctx *ctx,
        j = isl_basic_map_alloc_equality(ctx, bmap);
        if (j < 0)
                goto error;
+       isl_seq_clr(bmap->eq[j],
+                   1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
        isl_int_set_si(bmap->eq[j][1+bmap->nparam+input], -1);
        isl_int_set_si(bmap->eq[j][0], value);
+       bmap = isl_basic_map_simplify(ctx, bmap);
        return isl_basic_map_finalize(ctx, bmap);
 error:
        isl_basic_map_free(ctx, bmap);
@@ -1941,6 +2097,8 @@ error:
 struct isl_map *isl_basic_map_compute_divs(struct isl_ctx *ctx,
                struct isl_basic_map *bmap)
 {
+       if (!bmap)
+               return NULL;
        if (bmap->n_div == 0)
                return isl_map_from_basic_map(ctx, bmap);
        return isl_pip_basic_map_compute_divs(ctx, bmap);
@@ -2204,6 +2362,8 @@ struct isl_basic_set *isl_basic_map_deltas(struct isl_ctx *ctx,
                                            (struct isl_basic_map *)bset);
                if (j < 0)
                        goto error;
+               isl_seq_clr(bset->eq[j],
+                           1 + bset->nparam + bset->dim + bset->extra);
                isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
                isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
                isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
@@ -2350,6 +2510,8 @@ struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
                int j = isl_basic_map_alloc_equality(ctx, bmap);
                if (j < 0)
                        goto error;
+               isl_seq_clr(bmap->eq[j],
+                   1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
                isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
                isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
        }
@@ -2492,12 +2654,47 @@ int isl_basic_map_is_strict_subset(struct isl_ctx *ctx,
        return !is_subset;
 }
 
+static int basic_map_contains(struct isl_ctx *ctx, struct isl_basic_map *bmap,
+                                        struct isl_vec *vec)
+{
+       int i;
+       unsigned total;
+       isl_int s;
+
+       total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       if (total != vec->size)
+               return -1;
+
+       isl_int_init(s);
+
+       for (i = 0; i < bmap->n_eq; ++i) {
+               isl_seq_inner_product(vec->block.data, bmap->eq[i], total, &s);
+               if (!isl_int_is_zero(s)) {
+                       isl_int_clear(s);
+                       return 0;
+               }
+       }
+
+       for (i = 0; i < bmap->n_ineq; ++i) {
+               isl_seq_inner_product(vec->block.data, bmap->ineq[i], total, &s);
+               if (isl_int_is_neg(s)) {
+                       isl_int_clear(s);
+                       return 0;
+               }
+       }
+
+       isl_int_clear(s);
+
+       return 1;
+}
+
 int isl_basic_map_is_empty(struct isl_ctx *ctx,
                struct isl_basic_map *bmap)
 {
        struct isl_basic_set *bset = NULL;
-       struct isl_set *min = NULL;
+       struct isl_vec *sample = NULL;
        int empty;
+       unsigned total;
 
        if (!bmap)
                return -1;
@@ -2505,15 +2702,26 @@ int isl_basic_map_is_empty(struct isl_ctx *ctx,
        if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
                return 1;
 
-       bset = isl_basic_set_from_basic_map(ctx,
+       total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
+       if (bmap->sample && bmap->sample->size == total) {
+               int contains = basic_map_contains(ctx, bmap, bmap->sample);
+               if (contains < 0)
+                       return -1;
+               if (contains)
+                       return 0;
+       }
+       bset = isl_basic_map_underlying_set(ctx,
                        isl_basic_map_copy(ctx, bmap));
        if (!bset)
                return -1;
-       min = isl_basic_set_lexmin(ctx, bset);
-       if (!min)
+       sample = isl_basic_set_sample(ctx, bset);
+       if (!sample)
                return -1;
-       empty = min->n == 0;
-       isl_set_free(ctx, min);
+       empty = sample->size == 0;
+       if (bmap->sample)
+               isl_vec_free(ctx, bmap->sample);
+       bmap->sample = sample;
+
        return empty;
 }
 
@@ -2599,7 +2807,7 @@ static int add_div_constraints(struct isl_ctx *ctx,
 {
        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;
+       unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
 
        i = isl_basic_map_alloc_inequality(ctx, bmap);
        if (i < 0)
@@ -2637,6 +2845,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);
+                       isl_seq_clr(dst->div[j]+1+1+total,
+                                                   dst->extra - src->n_div);
                        if (add_div_constraints(ctx, dst, j) < 0)
                                goto error;
                }
@@ -2649,14 +2859,32 @@ error:
        return NULL;
 }
 
+struct isl_map *isl_map_align_divs(struct isl_ctx *ctx, struct isl_map *map)
+{
+       int i;
+
+       map = isl_map_compute_divs(ctx, map);
+       map = isl_map_cow(ctx, map);
+       if (!map)
+               return NULL;
+
+       for (i = 1; i < map->n; ++i)
+               map->p[0] = isl_basic_map_align_divs(ctx, map->p[0], map->p[i]);
+       for (i = 1; i < map->n; ++i)
+               map->p[i] = isl_basic_map_align_divs(ctx, map->p[i], map->p[0]);
+
+       return map;
+}
+
 static struct isl_map *add_cut_constraint(struct isl_ctx *ctx,
                struct isl_map *dst,
                struct isl_basic_map *src, isl_int *c,
-               unsigned len, unsigned extra, int oppose)
+               unsigned len, int oppose)
 {
        struct isl_basic_map *copy = NULL;
        int is_empty;
        int k;
+       unsigned total;
 
        copy = isl_basic_map_copy(ctx, src);
        copy = isl_basic_map_cow(ctx, copy);
@@ -2671,7 +2899,8 @@ static struct isl_map *add_cut_constraint(struct isl_ctx *ctx,
                isl_seq_neg(copy->ineq[k], c, len);
        else
                isl_seq_cpy(copy->ineq[k], c, len);
-       isl_seq_clr(copy->ineq[k]+len, extra);
+       total = 1 + copy->nparam + copy->n_in + copy->n_out + copy->extra;
+       isl_seq_clr(copy->ineq[k]+len, total - len);
        isl_inequality_negate(ctx, copy, k);
        copy = isl_basic_map_simplify(ctx, copy);
        copy = isl_basic_map_finalize(ctx, copy);
@@ -2724,14 +2953,12 @@ static struct isl_map *subtract(struct isl_ctx *ctx, struct isl_map *map,
        for (i = 0; i < bmap->n_eq; ++i) {
                for (j = 0; j < map->n; ++j) {
                        rest = add_cut_constraint(ctx, rest,
-                               map->p[j], bmap->eq[i],
-                               1+total, map->p[j]->n_div - bmap->n_div, 0);
+                               map->p[j], bmap->eq[i], 1+total, 0);
                        if (!rest)
                                goto error;
 
                        rest = add_cut_constraint(ctx, rest,
-                               map->p[j], bmap->eq[i],
-                               1+total, map->p[j]->n_div - bmap->n_div, 1);
+                               map->p[j], bmap->eq[i], 1+total, 1);
                        if (!rest)
                                goto error;
 
@@ -2745,15 +2972,14 @@ static struct isl_map *subtract(struct isl_ctx *ctx, struct isl_map *map,
                                goto error;
                        isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
                        isl_seq_clr(map->p[j]->eq[k]+1+total,
-                                       map->p[j]->n_div - bmap->n_div);
+                                       map->p[j]->extra - bmap->n_div);
                }
        }
 
        for (i = 0; i < bmap->n_ineq; ++i) {
                for (j = 0; j < map->n; ++j) {
                        rest = add_cut_constraint(ctx, rest,
-                               map->p[j], bmap->ineq[i],
-                               1+total, map->p[j]->n_div - bmap->n_div, 0);
+                               map->p[j], bmap->ineq[i], 1+total, 0);
                        if (!rest)
                                goto error;
 
@@ -2767,7 +2993,7 @@ static struct isl_map *subtract(struct isl_ctx *ctx, struct isl_map *map,
                                goto error;
                        isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
                        isl_seq_clr(map->p[j]->ineq[k]+1+total,
-                                       map->p[j]->n_div - bmap->n_div);
+                                       map->p[j]->extra - bmap->n_div);
                }
        }
 
@@ -2831,3 +3057,33 @@ error:
        isl_map_free(ctx, map);
        return NULL;
 }
+
+/* There is no need to cow as removing empty parts doesn't change
+ * the meaning of the set.
+ */
+struct isl_map *isl_map_remove_empty_parts(struct isl_ctx *ctx,
+       struct isl_map *map)
+{
+       int i;
+
+       if (!map)
+               return NULL;
+
+       for (i = map->n-1; i >= 0; --i) {
+               if (!F_ISSET(map->p[i], ISL_BASIC_MAP_EMPTY))
+                       continue;
+               isl_basic_map_free(ctx, map->p[i]);
+               if (i != map->n-1)
+                       map->p[i] = map->p[map->n-1];
+               map->n--;
+       }
+
+       return map;
+}
+
+struct isl_set *isl_set_remove_empty_parts(struct isl_ctx *ctx,
+       struct isl_set *set)
+{
+       return (struct isl_set *)
+               isl_map_remove_empty_parts(ctx, (struct isl_map *)set);
+}