isl_map_subtract.c

   1 #include "isl_seq.h"
   2 #include "isl_set.h"
   3 #include "isl_map.h"
   4 #include "isl_map_private.h"
   5
   6 static struct isl_map *add_cut_constraint(struct isl_map *dst,
   7                 struct isl_basic_map *src, isl_int *c,
   8                 unsigned len, int oppose)
   9 {
  10         struct isl_basic_map *copy = NULL;
  11         int is_empty;
  12         int k;
  13         unsigned total;
  14
  15         copy = isl_basic_map_copy(src);
  16         copy = isl_basic_map_cow(copy);
  17         if (!copy)
  18                 goto error;
  19         copy = isl_basic_map_extend_constraints(copy, 0, 1);
  20         k = isl_basic_map_alloc_inequality(copy);
  21         if (k < 0)
  22                 goto error;
  23         if (oppose)
  24                 isl_seq_neg(copy->ineq[k], c, len);
  25         else
  26                 isl_seq_cpy(copy->ineq[k], c, len);
  27         total = 1 + isl_basic_map_total_dim(copy);
  28         isl_seq_clr(copy->ineq[k]+len, total - len);
  29         isl_inequality_negate(copy, k);
  30         copy = isl_basic_map_simplify(copy);
  31         copy = isl_basic_map_finalize(copy);
  32         is_empty = isl_basic_map_is_empty(copy);
  33         if (is_empty < 0)
  34                 goto error;
  35         if (!is_empty)
  36                 dst = isl_map_add(dst, copy);
  37         else
  38                 isl_basic_map_free(copy);
  39         return dst;
  40 error:
  41         isl_basic_map_free(copy);
  42         isl_map_free(dst);
  43         return NULL;
  44 }
  45
  46 static struct isl_map *subtract(struct isl_map *map, struct isl_basic_map *bmap)
  47 {
  48         int i, j, k;
  49         unsigned flags = 0;
  50         struct isl_map *rest = NULL;
  51         unsigned max;
  52         unsigned total = isl_basic_map_total_dim(bmap);
  53
  54         map = isl_map_cow(map);
  55
  56         if (!map || !bmap)
  57                 goto error;
  58
  59         if (ISL_F_ISSET(map, ISL_MAP_DISJOINT))
  60                 ISL_FL_SET(flags, ISL_MAP_DISJOINT);
  61
  62         max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
  63         rest = isl_map_alloc_dim(isl_dim_copy(map->dim), max, flags);
  64         if (!rest)
  65                 goto error;
  66
  67         for (i = 0; i < map->n; ++i) {
  68                 map->p[i] = isl_basic_map_align_divs(map->p[i], bmap);
  69                 if (!map->p[i])
  70                         goto error;
  71         }
  72
  73         for (j = 0; j < map->n; ++j)
  74                 map->p[j] = isl_basic_map_cow(map->p[j]);
  75
  76         for (i = 0; i < bmap->n_eq; ++i) {
  77                 for (j = 0; j < map->n; ++j) {
  78                         rest = add_cut_constraint(rest,
  79                                 map->p[j], bmap->eq[i], 1+total, 0);
  80                         if (!rest)
  81                                 goto error;
  82
  83                         rest = add_cut_constraint(rest,
  84                                 map->p[j], bmap->eq[i], 1+total, 1);
  85                         if (!rest)
  86                                 goto error;
  87
  88                         map->p[j] = isl_basic_map_extend_constraints(map->p[j],
  89                                 1, 0);
  90                         if (!map->p[j])
  91                                 goto error;
  92                         k = isl_basic_map_alloc_equality(map->p[j]);
  93                         if (k < 0)
  94                                 goto error;
  95                         isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
  96                         isl_seq_clr(map->p[j]->eq[k]+1+total,
  97                                         map->p[j]->n_div - bmap->n_div);
  98                 }
  99         }
 100
 101         for (i = 0; i < bmap->n_ineq; ++i) {
 102                 for (j = 0; j < map->n; ++j) {
 103                         rest = add_cut_constraint(rest,
 104                                 map->p[j], bmap->ineq[i], 1+total, 0);
 105                         if (!rest)
 106                                 goto error;
 107
 108                         map->p[j] = isl_basic_map_extend_constraints(map->p[j],
 109                                 0, 1);
 110                         if (!map->p[j])
 111                                 goto error;
 112                         k = isl_basic_map_alloc_inequality(map->p[j]);
 113                         if (k < 0)
 114                                 goto error;
 115                         isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
 116                         isl_seq_clr(map->p[j]->ineq[k]+1+total,
 117                                         map->p[j]->n_div - bmap->n_div);
 118                 }
 119         }
 120
 121         isl_map_free(map);
 122         return rest;
 123 error:
 124         isl_map_free(map);
 125         isl_map_free(rest);
 126         return NULL;
 127 }
 128
 129 struct isl_map *isl_map_subtract(struct isl_map *map1, struct isl_map *map2)
 130 {
 131         int i;
 132         if (!map1 || !map2)
 133                 goto error;
 134
 135         isl_assert(map1->ctx, isl_dim_equal(map1->dim, map2->dim), goto error);
 136
 137         if (isl_map_is_empty(map2)) {
 138                 isl_map_free(map2);
 139                 return map1;
 140         }
 141
 142         map1 = isl_map_compute_divs(map1);
 143         map2 = isl_map_compute_divs(map2);
 144         if (!map1 || !map2)
 145                 goto error;
 146
 147         map1 = isl_map_remove_empty_parts(map1);
 148         map2 = isl_map_remove_empty_parts(map2);
 149
 150         for (i = 0; map1 && i < map2->n; ++i)
 151                 map1 = subtract(map1, map2->p[i]);
 152
 153         isl_map_free(map2);
 154         return map1;
 155 error:
 156         isl_map_free(map1);
 157         isl_map_free(map2);
 158         return NULL;
 159 }
 160
 161 struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
 162 {
 163         return (struct isl_set *)
 164                 isl_map_subtract(
 165                         (struct isl_map *)set1, (struct isl_map *)set2);
 166 }
 167
 168 /* Return 1 if "bmap" contains a single element.
 169  */
 170 int isl_basic_map_is_singleton(__isl_keep isl_basic_map *bmap)
 171 {
 172         if (!bmap)
 173                 return -1;
 174         if (bmap->n_div)
 175                 return 0;
 176         if (bmap->n_ineq)
 177                 return 0;
 178         return bmap->n_eq == isl_basic_map_total_dim(bmap);
 179 }
 180
 181 /* Return 1 if "map" contains a single element.
 182  */
 183 int isl_map_is_singleton(__isl_keep isl_map *map)
 184 {
 185         if (!map)
 186                 return -1;
 187         if (map->n != 1)
 188                 return 0;
 189
 190         return isl_basic_map_is_singleton(map->p[0]);
 191 }
 192
 193 /* Given a singleton basic map, extract the single element
 194  * as an isl_vec.
 195  */
 196 static __isl_give isl_vec *singleton_extract_point(__isl_keep isl_basic_map *bmap)
 197 {
 198         int i, j;
 199         unsigned dim;
 200         struct isl_vec *point;
 201         isl_int m;
 202
 203         if (!bmap)
 204                 return NULL;
 205
 206         dim = isl_basic_map_total_dim(bmap);
 207         isl_assert(bmap->ctx, bmap->n_eq == dim, return NULL);
 208         point = isl_vec_alloc(bmap->ctx, 1 + dim);
 209         if (!point)
 210                 return NULL;
 211
 212         isl_int_init(m);
 213
 214         isl_int_set_si(point->el[0], 1);
 215         for (j = 0; j < bmap->n_eq; ++j) {
 216                 int s;
 217                 int i = dim - 1 - j;
 218                 isl_assert(bmap->ctx,
 219                     isl_seq_first_non_zero(bmap->eq[j] + 1, i) == -1,
 220                     goto error);
 221                 isl_assert(bmap->ctx,
 222                     isl_int_is_one(bmap->eq[j][1 + i]) ||
 223                     isl_int_is_negone(bmap->eq[j][1 + i]),
 224                     goto error);
 225                 isl_assert(bmap->ctx,
 226                     isl_seq_first_non_zero(bmap->eq[j]+1+i+1, dim-i-1) == -1,
 227                     goto error);
 228
 229                 isl_int_gcd(m, point->el[0], bmap->eq[j][1 + i]);
 230                 isl_int_divexact(m, bmap->eq[j][1 + i], m);
 231                 isl_int_abs(m, m);
 232                 isl_seq_scale(point->el, point->el, m, 1 + i);
 233                 isl_int_divexact(m, point->el[0], bmap->eq[j][1 + i]);
 234                 isl_int_neg(m, m);
 235                 isl_int_mul(point->el[1 + i], m, bmap->eq[j][0]);
 236         }
 237
 238         isl_int_clear(m);
 239         return point;
 240 error:
 241         isl_int_clear(m);
 242         isl_vec_free(point);
 243         return NULL;
 244 }
 245
 246 /* Return 1 if "bmap" contains the point "point".
 247  * "bmap" is assumed to have known divs.
 248  * The point is first extended with the divs and then passed
 249  * to basic_map_contains.
 250  */
 251 static int basic_map_contains_point(__isl_keep isl_basic_map *bmap,
 252         __isl_keep isl_vec *point)
 253 {
 254         int i;
 255         struct isl_vec *vec;
 256         unsigned dim;
 257         int contains;
 258
 259         if (!bmap || !point)
 260                 return -1;
 261         if (bmap->n_div == 0)
 262                 return isl_basic_map_contains(bmap, point);
 263
 264         dim = isl_basic_map_total_dim(bmap) - bmap->n_div;
 265         vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
 266         if (!vec)
 267                 return -1;
 268
 269         isl_seq_cpy(vec->el, point->el, point->size);
 270         for (i = 0; i < bmap->n_div; ++i) {
 271                 isl_seq_inner_product(bmap->div[i] + 1, vec->el,
 272                                         1 + dim + i, &vec->el[1+dim+i]);
 273                 isl_int_fdiv_q(vec->el[1+dim+i], vec->el[1+dim+i],
 274                                 bmap->div[i][0]);
 275         }
 276
 277         contains = isl_basic_map_contains(bmap, vec);
 278
 279         isl_vec_free(vec);
 280         return contains;
 281 }
 282
 283 /* Return 1 is the singleton map "map1" is a subset of "map2",
 284  * i.e., if the single element of "map1" is also an element of "map2".
 285  */
 286 static int map_is_singleton_subset(__isl_keep isl_map *map1,
 287         __isl_keep isl_map *map2)
 288 {
 289         int i;
 290         int is_subset = 0;
 291         struct isl_vec *point;
 292
 293         if (!map1 || !map2)
 294                 return -1;
 295         if (map1->n != 1)
 296                 return -1;
 297
 298         point = singleton_extract_point(map1->p[0]);
 299         if (!point)
 300                 return -1;
 301
 302         for (i = 0; i < map2->n; ++i) {
 303                 is_subset = basic_map_contains_point(map2->p[i], point);
 304                 if (is_subset)
 305                         break;
 306         }
 307
 308         isl_vec_free(point);
 309         return is_subset;
 310 }
 311
 312 int isl_map_is_subset(struct isl_map *map1, struct isl_map *map2)
 313 {
 314         int is_subset = 0;
 315         struct isl_map *diff;
 316
 317         if (!map1 || !map2)
 318                 return -1;
 319
 320         if (isl_map_is_empty(map1))
 321                 return 1;
 322
 323         if (isl_map_is_empty(map2))
 324                 return 0;
 325
 326         if (isl_map_fast_is_universe(map2))
 327                 return 1;
 328
 329         map2 = isl_map_compute_divs(isl_map_copy(map2));
 330         if (isl_map_is_singleton(map1)) {
 331                 is_subset = map_is_singleton_subset(map1, map2);
 332                 isl_map_free(map2);
 333                 return is_subset;
 334         }
 335         diff = isl_map_subtract(isl_map_copy(map1), map2);
 336         if (!diff)
 337                 return -1;
 338
 339         is_subset = isl_map_is_empty(diff);
 340         isl_map_free(diff);
 341
 342         return is_subset;
 343 }
 344
 345 int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
 346 {
 347         return isl_map_is_subset(
 348                         (struct isl_map *)set1, (struct isl_map *)set2);
 349 }