isl_map.c

   1 #include <string.h>
   2 #include <strings.h>
   3 #include "isl_ctx.h"
   4 #include "isl_blk.h"
   5 #include "isl_lp.h"
   6 #include "isl_seq.h"
   7 #include "isl_set.h"
   8 #include "isl_map.h"
   9 #include "isl_map_private.h"
  10 #include "isl_map_piplib.h"
  11 #include "isl_sample.h"
  12 #include "isl_vec.h"
  13
  14 static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
  15                 struct isl_basic_map *bmap,
  16                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
  17                 unsigned n_eq, unsigned n_ineq)
  18 {
  19         int i;
  20         size_t row_size = 1 + nparam + n_in + n_out + extra;
  21
  22         bmap->block = isl_blk_alloc(ctx, (n_eq + n_ineq) * row_size);
  23         if (isl_blk_is_error(bmap->block)) {
  24                 free(bmap);
  25                 return NULL;
  26         }
  27
  28         bmap->eq = isl_alloc_array(ctx, isl_int *, n_eq + n_ineq);
  29         if (!bmap->eq) {
  30                 isl_blk_free(ctx, bmap->block);
  31                 free(bmap);
  32                 return NULL;
  33         }
  34
  35         if (extra == 0) {
  36                 bmap->block2 = isl_blk_empty();
  37                 bmap->div = NULL;
  38         } else {
  39                 bmap->block2 = isl_blk_alloc(ctx, extra * (1 + row_size));
  40                 if (isl_blk_is_error(bmap->block2)) {
  41                         free(bmap->eq);
  42                         isl_blk_free(ctx, bmap->block);
  43                         free(bmap);
  44                         return NULL;
  45                 }
  46
  47                 bmap->div = isl_alloc_array(ctx, isl_int *, extra);
  48                 if (!bmap->div) {
  49                         isl_blk_free(ctx, bmap->block2);
  50                         free(bmap->eq);
  51                         isl_blk_free(ctx, bmap->block);
  52                         free(bmap);
  53                         return NULL;
  54                 }
  55         }
  56
  57         for (i = 0; i < n_eq + n_ineq; ++i)
  58                 bmap->eq[i] = bmap->block.data + i * row_size;
  59
  60         for (i = 0; i < extra; ++i)
  61                 bmap->div[i] = bmap->block2.data + i * (1 + row_size);
  62
  63         bmap->ctx = ctx;
  64         isl_ctx_ref(ctx);
  65         bmap->ref = 1;
  66         bmap->flags = 0;
  67         bmap->c_size = n_eq + n_ineq;
  68         bmap->ineq = bmap->eq + n_eq;
  69         bmap->nparam = nparam;
  70         bmap->n_in = n_in;
  71         bmap->n_out = n_out;
  72         bmap->extra = extra;
  73         bmap->n_eq = 0;
  74         bmap->n_ineq = 0;
  75         bmap->n_div = 0;
  76         bmap->sample = NULL;
  77
  78         return bmap;
  79 }
  80
  81 struct isl_basic_set *isl_basic_set_alloc(struct isl_ctx *ctx,
  82                 unsigned nparam, unsigned dim, unsigned extra,
  83                 unsigned n_eq, unsigned n_ineq)
  84 {
  85         struct isl_basic_map *bmap;
  86         bmap = isl_basic_map_alloc(ctx, nparam, 0, dim, extra, n_eq, n_ineq);
  87         return (struct isl_basic_set *)bmap;
  88 }
  89
  90 struct isl_basic_map *isl_basic_map_alloc(struct isl_ctx *ctx,
  91                 unsigned nparam, unsigned in, unsigned out, unsigned extra,
  92                 unsigned n_eq, unsigned n_ineq)
  93 {
  94         struct isl_basic_map *bmap;
  95
  96         bmap = isl_alloc_type(ctx, struct isl_basic_map);
  97         if (!bmap)
  98                 return NULL;
  99
 100         return basic_map_init(ctx, bmap,
 101                                nparam, in, out, extra, n_eq, n_ineq);
 102 }
 103
 104 static void dup_constraints(
 105                 struct isl_basic_map *dst, struct isl_basic_map *src)
 106 {
 107         int i;
 108         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
 109
 110         for (i = 0; i < src->n_eq; ++i) {
 111                 int j = isl_basic_map_alloc_equality(dst);
 112                 isl_seq_cpy(dst->eq[j], src->eq[i], 1+total);
 113         }
 114
 115         for (i = 0; i < src->n_ineq; ++i) {
 116                 int j = isl_basic_map_alloc_inequality(dst);
 117                 isl_seq_cpy(dst->ineq[j], src->ineq[i], 1+total);
 118         }
 119
 120         for (i = 0; i < src->n_div; ++i) {
 121                 int j = isl_basic_map_alloc_div(dst);
 122                 isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
 123         }
 124         F_SET(dst, ISL_BASIC_SET_FINAL);
 125 }
 126
 127 struct isl_basic_map *isl_basic_map_dup(struct isl_basic_map *bmap)
 128 {
 129         struct isl_basic_map *dup;
 130
 131         if (!bmap)
 132                 return NULL;
 133         dup = isl_basic_map_alloc(bmap->ctx, bmap->nparam,
 134                         bmap->n_in, bmap->n_out,
 135                         bmap->n_div, bmap->n_eq, bmap->n_ineq);
 136         if (!dup)
 137                 return NULL;
 138         dup->flags = bmap->flags;
 139         dup_constraints(dup, bmap);
 140         dup->sample = isl_vec_copy(bmap->ctx, bmap->sample);
 141         return dup;
 142 }
 143
 144 struct isl_basic_set *isl_basic_set_dup(struct isl_basic_set *bset)
 145 {
 146         struct isl_basic_map *dup;
 147
 148         dup = isl_basic_map_dup((struct isl_basic_map *)bset);
 149         return (struct isl_basic_set *)dup;
 150 }
 151
 152 struct isl_basic_set *isl_basic_set_copy(struct isl_basic_set *bset)
 153 {
 154         if (!bset)
 155                 return NULL;
 156
 157         if (F_ISSET(bset, ISL_BASIC_SET_FINAL)) {
 158                 bset->ref++;
 159                 return bset;
 160         }
 161         return isl_basic_set_dup(bset);
 162 }
 163
 164 struct isl_set *isl_set_copy(struct isl_set *set)
 165 {
 166         if (!set)
 167                 return NULL;
 168
 169         set->ref++;
 170         return set;
 171 }
 172
 173 struct isl_basic_map *isl_basic_map_copy(struct isl_basic_map *bmap)
 174 {
 175         if (!bmap)
 176                 return NULL;
 177
 178         if (F_ISSET(bmap, ISL_BASIC_SET_FINAL)) {
 179                 bmap->ref++;
 180                 return bmap;
 181         }
 182         return isl_basic_map_dup(bmap);
 183 }
 184
 185 struct isl_map *isl_map_copy(struct isl_map *map)
 186 {
 187         if (!map)
 188                 return NULL;
 189
 190         map->ref++;
 191         return map;
 192 }
 193
 194 void isl_basic_map_free(struct isl_basic_map *bmap)
 195 {
 196         if (!bmap)
 197                 return;
 198
 199         if (--bmap->ref > 0)
 200                 return;
 201
 202         isl_ctx_deref(bmap->ctx);
 203         free(bmap->div);
 204         isl_blk_free(bmap->ctx, bmap->block2);
 205         free(bmap->eq);
 206         isl_blk_free(bmap->ctx, bmap->block);
 207         isl_vec_free(bmap->ctx, bmap->sample);
 208         free(bmap);
 209 }
 210
 211 void isl_basic_set_free(struct isl_basic_set *bset)
 212 {
 213         isl_basic_map_free((struct isl_basic_map *)bset);
 214 }
 215
 216 int isl_basic_map_alloc_equality(struct isl_basic_map *bmap)
 217 {
 218         struct isl_ctx *ctx;
 219         if (!bmap)
 220                 return -1;
 221         ctx = bmap->ctx;
 222         isl_assert(ctx, bmap->n_eq + bmap->n_ineq < bmap->c_size, return -1);
 223         isl_assert(ctx, bmap->eq + bmap->n_eq <= bmap->ineq, return -1);
 224         if (bmap->eq + bmap->n_eq == bmap->ineq) {
 225                 isl_int *t;
 226                 int j = isl_basic_map_alloc_inequality(bmap);
 227                 if (j < 0)
 228                         return -1;
 229                 t = bmap->ineq[0];
 230                 bmap->ineq[0] = bmap->ineq[j];
 231                 bmap->ineq[j] = t;
 232                 bmap->n_ineq--;
 233                 bmap->ineq++;
 234         }
 235         return bmap->n_eq++;
 236 }
 237
 238 int isl_basic_set_alloc_equality(struct isl_basic_set *bset)
 239 {
 240         return isl_basic_map_alloc_equality((struct isl_basic_map *)bset);
 241 }
 242
 243 int isl_basic_map_free_equality(struct isl_basic_map *bmap, unsigned n)
 244 {
 245         if (!bmap)
 246                 return -1;
 247         isl_assert(bmap->ctx, n <= bmap->n_eq, return -1);
 248         bmap->n_eq -= n;
 249         return 0;
 250 }
 251
 252 int isl_basic_map_drop_equality(struct isl_basic_map *bmap, unsigned pos)
 253 {
 254         isl_int *t;
 255         if (!bmap)
 256                 return -1;
 257         isl_assert(bmap->ctx, pos < bmap->n_eq, return -1);
 258
 259         if (pos != bmap->n_eq - 1) {
 260                 t = bmap->eq[pos];
 261                 bmap->eq[pos] = bmap->eq[bmap->n_eq - 1];
 262                 bmap->eq[bmap->n_eq - 1] = t;
 263         }
 264         bmap->n_eq--;
 265         return 0;
 266 }
 267
 268 int isl_basic_set_drop_equality(struct isl_basic_set *bset, unsigned pos)
 269 {
 270         return isl_basic_map_drop_equality((struct isl_basic_map *)bset, pos);
 271 }
 272
 273 void isl_basic_map_inequality_to_equality(
 274                 struct isl_basic_map *bmap, unsigned pos)
 275 {
 276         isl_int *t;
 277
 278         t = bmap->ineq[pos];
 279         bmap->ineq[pos] = bmap->ineq[0];
 280         bmap->ineq[0] = bmap->eq[bmap->n_eq];
 281         bmap->eq[bmap->n_eq] = t;
 282         bmap->n_eq++;
 283         bmap->n_ineq--;
 284         bmap->ineq++;
 285 }
 286
 287 int isl_basic_map_alloc_inequality(struct isl_basic_map *bmap)
 288 {
 289         struct isl_ctx *ctx;
 290         if (!bmap)
 291                 return -1;
 292         ctx = bmap->ctx;
 293         isl_assert(ctx, (bmap->ineq - bmap->eq) + bmap->n_ineq < bmap->c_size,
 294                         return -1);
 295         F_CLR(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
 296         F_CLR(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
 297         return bmap->n_ineq++;
 298 }
 299
 300 int isl_basic_set_alloc_inequality(struct isl_basic_set *bset)
 301 {
 302         return isl_basic_map_alloc_inequality((struct isl_basic_map *)bset);
 303 }
 304
 305 int isl_basic_map_free_inequality(struct isl_basic_map *bmap, unsigned n)
 306 {
 307         if (!bmap)
 308                 return -1;
 309         isl_assert(bmap->ctx, n <= bmap->n_ineq, return -1);
 310         bmap->n_ineq -= n;
 311         return 0;
 312 }
 313
 314 int isl_basic_set_free_inequality(struct isl_basic_set *bset, unsigned n)
 315 {
 316         return isl_basic_map_free_inequality((struct isl_basic_map *)bset, n);
 317 }
 318
 319 int isl_basic_map_drop_inequality(struct isl_basic_map *bmap, unsigned pos)
 320 {
 321         isl_int *t;
 322         if (!bmap)
 323                 return -1;
 324         isl_assert(bmap->ctx, pos < bmap->n_ineq, return -1);
 325
 326         if (pos != bmap->n_ineq - 1) {
 327                 t = bmap->ineq[pos];
 328                 bmap->ineq[pos] = bmap->ineq[bmap->n_ineq - 1];
 329                 bmap->ineq[bmap->n_ineq - 1] = t;
 330         }
 331         bmap->n_ineq--;
 332         return 0;
 333 }
 334
 335 int isl_basic_set_drop_inequality(struct isl_basic_set *bset, unsigned pos)
 336 {
 337         return isl_basic_map_drop_inequality((struct isl_basic_map *)bset, pos);
 338 }
 339
 340 int isl_basic_map_alloc_div(struct isl_basic_map *bmap)
 341 {
 342         if (!bmap)
 343                 return -1;
 344         isl_assert(bmap->ctx, bmap->n_div < bmap->extra, return -1);
 345         return bmap->n_div++;
 346 }
 347
 348 int isl_basic_map_free_div(struct isl_basic_map *bmap, unsigned n)
 349 {
 350         if (!bmap)
 351                 return -1;
 352         isl_assert(bmap->ctx, n <= bmap->n_div, return -1);
 353         bmap->n_div -= n;
 354         return 0;
 355 }
 356
 357 /* Copy constraint from src to dst, putting the vars of src at offset
 358  * dim_off in dst and the divs of src at offset div_off in dst.
 359  * If both sets are actually map, then dim_off applies to the input
 360  * variables.
 361  */
 362 static void copy_constraint(struct isl_basic_map *dst_map, isl_int *dst,
 363                             struct isl_basic_map *src_map, isl_int *src,
 364                             unsigned in_off, unsigned out_off, unsigned div_off)
 365 {
 366         isl_int_set(dst[0], src[0]);
 367         isl_seq_cpy(dst+1, src+1, isl_min(dst_map->nparam, src_map->nparam));
 368         if (dst_map->nparam > src_map->nparam)
 369                 isl_seq_clr(dst+1+src_map->nparam,
 370                                 dst_map->nparam - src_map->nparam);
 371         isl_seq_clr(dst+1+dst_map->nparam, in_off);
 372         isl_seq_cpy(dst+1+dst_map->nparam+in_off,
 373                     src+1+src_map->nparam,
 374                     isl_min(dst_map->n_in-in_off, src_map->n_in));
 375         if (dst_map->n_in-in_off > src_map->n_in)
 376                 isl_seq_clr(dst+1+dst_map->nparam+in_off+src_map->n_in,
 377                                 dst_map->n_in - in_off - src_map->n_in);
 378         isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in, out_off);
 379         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+out_off,
 380                     src+1+src_map->nparam+src_map->n_in,
 381                     isl_min(dst_map->n_out-out_off, src_map->n_out));
 382         if (dst_map->n_out-out_off > src_map->n_out)
 383                 isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+out_off+src_map->n_out,
 384                                 dst_map->n_out - out_off - src_map->n_out);
 385         isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out, div_off);
 386         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+div_off,
 387                     src+1+src_map->nparam+src_map->n_in+src_map->n_out,
 388                     isl_min(dst_map->extra-div_off, src_map->extra));
 389         if (dst_map->extra-div_off > src_map->extra)
 390                 isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+
 391                                 div_off+src_map->extra,
 392                                 dst_map->extra - div_off - src_map->extra);
 393 }
 394
 395 static void copy_div(struct isl_basic_map *dst_map, isl_int *dst,
 396                      struct isl_basic_map *src_map, isl_int *src,
 397                      unsigned in_off, unsigned out_off, unsigned div_off)
 398 {
 399         isl_int_set(dst[0], src[0]);
 400         copy_constraint(dst_map, dst+1, src_map, src+1, in_off, out_off, div_off);
 401 }
 402
 403 static struct isl_basic_map *add_constraints(struct isl_basic_map *bmap1,
 404                 struct isl_basic_map *bmap2, unsigned i_pos, unsigned o_pos)
 405 {
 406         int i;
 407         unsigned div_off = bmap1->n_div;
 408
 409         for (i = 0; i < bmap2->n_eq; ++i) {
 410                 int i1 = isl_basic_map_alloc_equality(bmap1);
 411                 if (i1 < 0)
 412                         goto error;
 413                 copy_constraint(bmap1, bmap1->eq[i1], bmap2, bmap2->eq[i],
 414                                 i_pos, o_pos, div_off);
 415         }
 416
 417         for (i = 0; i < bmap2->n_ineq; ++i) {
 418                 int i1 = isl_basic_map_alloc_inequality(bmap1);
 419                 if (i1 < 0)
 420                         goto error;
 421                 copy_constraint(bmap1, bmap1->ineq[i1], bmap2, bmap2->ineq[i],
 422                                 i_pos, o_pos, div_off);
 423         }
 424
 425         for (i = 0; i < bmap2->n_div; ++i) {
 426                 int i1 = isl_basic_map_alloc_div(bmap1);
 427                 if (i1 < 0)
 428                         goto error;
 429                 copy_div(bmap1, bmap1->div[i1], bmap2, bmap2->div[i],
 430                          i_pos, o_pos, div_off);
 431         }
 432
 433         isl_basic_map_free(bmap2);
 434
 435         return bmap1;
 436
 437 error:
 438         isl_basic_map_free(bmap1);
 439         isl_basic_map_free(bmap2);
 440         return NULL;
 441 }
 442
 443 struct isl_basic_map *isl_basic_map_extend(struct isl_basic_map *base,
 444                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
 445                 unsigned n_eq, unsigned n_ineq)
 446 {
 447         struct isl_basic_map *ext;
 448         unsigned flags;
 449         int dims_ok;
 450
 451         base = isl_basic_map_cow(base);
 452         if (!base)
 453                 goto error;
 454
 455         dims_ok = base && base->nparam == nparam &&
 456                   base->n_in == n_in && base->n_out == n_out &&
 457                   base->extra >= base->n_div + extra;
 458
 459         if (dims_ok && n_eq == 0 && n_ineq == 0)
 460                 return base;
 461
 462         isl_assert(base->ctx, base->nparam <= nparam, goto error);
 463         isl_assert(base->ctx, base->n_in <= n_in, goto error);
 464         isl_assert(base->ctx, base->n_out <= n_out, goto error);
 465         extra += base->extra;
 466         n_eq += base->n_eq;
 467         n_ineq += base->n_ineq;
 468
 469         ext = isl_basic_map_alloc(base->ctx, nparam, n_in, n_out,
 470                                         extra, n_eq, n_ineq);
 471         if (!base)
 472                 return ext;
 473         if (!ext)
 474                 goto error;
 475
 476         flags = base->flags;
 477         ext = add_constraints(ext, base, 0, 0);
 478         if (ext) {
 479                 ext->flags = flags;
 480                 F_CLR(ext, ISL_BASIC_SET_FINAL);
 481         }
 482
 483         return ext;
 484
 485 error:
 486         isl_basic_map_free(base);
 487         return NULL;
 488 }
 489
 490 struct isl_basic_set *isl_basic_set_extend(struct isl_basic_set *base,
 491                 unsigned nparam, unsigned dim, unsigned extra,
 492                 unsigned n_eq, unsigned n_ineq)
 493 {
 494         return (struct isl_basic_set *)
 495                 isl_basic_map_extend((struct isl_basic_map *)base,
 496                                         nparam, 0, dim, extra, n_eq, n_ineq);
 497 }
 498
 499 struct isl_basic_set *isl_basic_set_cow(struct isl_basic_set *bset)
 500 {
 501         return (struct isl_basic_set *)
 502                 isl_basic_map_cow((struct isl_basic_map *)bset);
 503 }
 504
 505 struct isl_basic_map *isl_basic_map_cow(struct isl_basic_map *bmap)
 506 {
 507         if (!bmap)
 508                 return NULL;
 509
 510         if (bmap->ref > 1) {
 511                 bmap->ref--;
 512                 bmap = isl_basic_map_dup(bmap);
 513         }
 514         F_CLR(bmap, ISL_BASIC_SET_FINAL);
 515         return bmap;
 516 }
 517
 518 struct isl_set *isl_set_cow(struct isl_set *set)
 519 {
 520         if (!set)
 521                 return NULL;
 522
 523         if (set->ref == 1)
 524                 return set;
 525         set->ref--;
 526         return isl_set_dup(set);
 527 }
 528
 529 struct isl_map *isl_map_cow(struct isl_map *map)
 530 {
 531         if (!map)
 532                 return NULL;
 533
 534         if (map->ref == 1)
 535                 return map;
 536         map->ref--;
 537         return isl_map_dup(map);
 538 }
 539
 540 static void swap_vars(struct isl_blk blk, isl_int *a,
 541                         unsigned a_len, unsigned b_len)
 542 {
 543         isl_seq_cpy(blk.data, a+a_len, b_len);
 544         isl_seq_cpy(blk.data+b_len, a, a_len);
 545         isl_seq_cpy(a, blk.data, b_len+a_len);
 546 }
 547
 548 struct isl_basic_set *isl_basic_set_swap_vars(
 549                 struct isl_basic_set *bset, unsigned n)
 550 {
 551         int i;
 552         struct isl_blk blk;
 553
 554         if (!bset)
 555                 goto error;
 556
 557         isl_assert(bset->ctx, n <= bset->dim, goto error);
 558
 559         if (n == bset->dim)
 560                 return bset;
 561
 562         bset = isl_basic_set_cow(bset);
 563         if (!bset)
 564                 return NULL;
 565
 566         blk = isl_blk_alloc(bset->ctx, bset->dim);
 567         if (isl_blk_is_error(blk))
 568                 goto error;
 569
 570         for (i = 0; i < bset->n_eq; ++i)
 571                 swap_vars(blk,
 572                           bset->eq[i]+1+bset->nparam, n, bset->dim - n);
 573
 574         for (i = 0; i < bset->n_ineq; ++i)
 575                 swap_vars(blk,
 576                           bset->ineq[i]+1+bset->nparam, n, bset->dim - n);
 577
 578         for (i = 0; i < bset->n_div; ++i)
 579                 swap_vars(blk,
 580                           bset->div[i]+1+1+bset->nparam, n, bset->dim - n);
 581
 582         isl_blk_free(bset->ctx, blk);
 583
 584         return bset;
 585
 586 error:
 587         isl_basic_set_free(bset);
 588         return NULL;
 589 }
 590
 591 struct isl_set *isl_set_swap_vars(struct isl_set *set, unsigned n)
 592 {
 593         int i;
 594         set = isl_set_cow(set);
 595         if (!set)
 596                 return NULL;
 597
 598         for (i = 0; i < set->n; ++i) {
 599                 set->p[i] = isl_basic_set_swap_vars(set->p[i], n);
 600                 if (!set->p[i]) {
 601                         isl_set_free(set);
 602                         return NULL;
 603                 }
 604         }
 605         return set;
 606 }
 607
 608 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
 609 {
 610         isl_seq_cpy(c, c + n, rem);
 611         isl_seq_clr(c + rem, n);
 612 }
 613
 614 struct isl_basic_set *isl_basic_set_drop_dims(
 615                 struct isl_basic_set *bset, unsigned first, unsigned n)
 616 {
 617         int i;
 618
 619         if (!bset)
 620                 goto error;
 621
 622         isl_assert(bset->ctx, first + n <= bset->dim, goto error);
 623
 624         if (n == 0)
 625                 return bset;
 626
 627         bset = isl_basic_set_cow(bset);
 628         if (!bset)
 629                 return NULL;
 630
 631         for (i = 0; i < bset->n_eq; ++i)
 632                 constraint_drop_vars(bset->eq[i]+1+bset->nparam+first, n,
 633                                      (bset->dim-first-n)+bset->extra);
 634
 635         for (i = 0; i < bset->n_ineq; ++i)
 636                 constraint_drop_vars(bset->ineq[i]+1+bset->nparam+first, n,
 637                                      (bset->dim-first-n)+bset->extra);
 638
 639         for (i = 0; i < bset->n_div; ++i)
 640                 constraint_drop_vars(bset->div[i]+1+1+bset->nparam+first, n,
 641                                      (bset->dim-first-n)+bset->extra);
 642
 643         bset->dim -= n;
 644         bset->extra += n;
 645
 646         return isl_basic_set_finalize(bset);
 647 error:
 648         isl_basic_set_free(bset);
 649         return NULL;
 650 }
 651
 652 static struct isl_set *isl_set_drop_dims(
 653                 struct isl_set *set, unsigned first, unsigned n)
 654 {
 655         int i;
 656
 657         if (!set)
 658                 goto error;
 659
 660         isl_assert(set->ctx, first + n <= set->dim, goto error);
 661
 662         if (n == 0)
 663                 return set;
 664         set = isl_set_cow(set);
 665         if (!set)
 666                 goto error;
 667
 668         for (i = 0; i < set->n; ++i) {
 669                 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
 670                 if (!set->p[i])
 671                         goto error;
 672         }
 673         set->dim -= n;
 674
 675         return set;
 676 error:
 677         isl_set_free(set);
 678         return NULL;
 679 }
 680
 681 /*
 682  * We don't cow, as the div is assumed to be redundant.
 683  */
 684 static struct isl_basic_map *isl_basic_map_drop_div(
 685                 struct isl_basic_map *bmap, unsigned div)
 686 {
 687         int i;
 688         unsigned pos;
 689
 690         if (!bmap)
 691                 goto error;
 692
 693         pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
 694
 695         isl_assert(bmap->ctx, div < bmap->n_div, goto error);
 696
 697         for (i = 0; i < bmap->n_eq; ++i)
 698                 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
 699
 700         for (i = 0; i < bmap->n_ineq; ++i) {
 701                 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
 702                         isl_basic_map_drop_inequality(bmap, i);
 703                         --i;
 704                         continue;
 705                 }
 706                 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
 707         }
 708
 709         for (i = 0; i < bmap->n_div; ++i)
 710                 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
 711
 712         if (div != bmap->n_div - 1) {
 713                 int j;
 714                 isl_int *t = bmap->div[div];
 715
 716                 for (j = div; j < bmap->n_div - 1; ++j)
 717                         bmap->div[j] = bmap->div[j+1];
 718
 719                 bmap->div[bmap->n_div - 1] = t;
 720         }
 721         isl_basic_map_free_div(bmap, 1);
 722
 723         return bmap;
 724 error:
 725         isl_basic_map_free(bmap);
 726         return NULL;
 727 }
 728
 729 static void eliminate_div(struct isl_basic_map *bmap, isl_int *eq, unsigned div)
 730 {
 731         int i;
 732         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
 733         unsigned len;
 734         len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 735
 736         for (i = 0; i < bmap->n_eq; ++i)
 737                 if (bmap->eq[i] != eq)
 738                         isl_seq_elim(bmap->eq[i], eq, pos, len, NULL);
 739
 740         for (i = 0; i < bmap->n_ineq; ++i)
 741                 isl_seq_elim(bmap->ineq[i], eq, pos, len, NULL);
 742
 743         /* We need to be careful about circular definitions,
 744          * so for now we just remove the definitions of other divs that
 745          * depend on this div and (possibly) recompute them later.
 746          */
 747         for (i = 0; i < bmap->n_div; ++i)
 748                 if (!isl_int_is_zero(bmap->div[i][0]) &&
 749                     !isl_int_is_zero(bmap->div[i][1 + pos]))
 750                         isl_seq_clr(bmap->div[i], 1 + len);
 751
 752         isl_basic_map_drop_div(bmap, div);
 753 }
 754
 755 struct isl_basic_map *isl_basic_map_set_to_empty(struct isl_basic_map *bmap)
 756 {
 757         int i = 0;
 758         unsigned total;
 759         if (!bmap)
 760                 goto error;
 761         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
 762         isl_basic_map_free_div(bmap, bmap->n_div);
 763         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 764         if (bmap->n_eq > 0)
 765                 isl_basic_map_free_equality(bmap, bmap->n_eq-1);
 766         else {
 767                 isl_basic_map_alloc_equality(bmap);
 768                 if (i < 0)
 769                         goto error;
 770         }
 771         isl_int_set_si(bmap->eq[i][0], 1);
 772         isl_seq_clr(bmap->eq[i]+1, total);
 773         F_SET(bmap, ISL_BASIC_MAP_EMPTY);
 774         return isl_basic_map_finalize(bmap);
 775 error:
 776         isl_basic_map_free(bmap);
 777         return NULL;
 778 }
 779
 780 struct isl_basic_set *isl_basic_set_set_to_empty(struct isl_basic_set *bset)
 781 {
 782         return (struct isl_basic_set *)
 783                 isl_basic_map_set_to_empty((struct isl_basic_map *)bset);
 784 }
 785
 786 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
 787 {
 788         isl_int *t = bmap->eq[a];
 789         bmap->eq[a] = bmap->eq[b];
 790         bmap->eq[b] = t;
 791 }
 792
 793 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
 794 {
 795         isl_int *t = bmap->ineq[a];
 796         bmap->ineq[a] = bmap->ineq[b];
 797         bmap->ineq[b] = t;
 798 }
 799
 800 static void swap_div(struct isl_basic_map *bmap, int a, int b)
 801 {
 802         int i;
 803         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
 804         isl_int *t = bmap->div[a];
 805         bmap->div[a] = bmap->div[b];
 806         bmap->div[b] = t;
 807
 808         for (i = 0; i < bmap->n_eq; ++i)
 809                 isl_int_swap(bmap->eq[i][1+off+a], bmap->eq[i][1+off+b]);
 810
 811         for (i = 0; i < bmap->n_ineq; ++i)
 812                 isl_int_swap(bmap->ineq[i][1+off+a], bmap->ineq[i][1+off+b]);
 813
 814         for (i = 0; i < bmap->n_div; ++i)
 815                 isl_int_swap(bmap->div[i][1+1+off+a], bmap->div[i][1+1+off+b]);
 816 }
 817
 818 struct isl_basic_map *isl_basic_map_gauss(
 819         struct isl_basic_map *bmap, int *progress)
 820 {
 821         int k;
 822         int done;
 823         int last_var;
 824         unsigned total_var;
 825         unsigned total;
 826
 827         if (!bmap)
 828                 return NULL;
 829
 830         total_var = bmap->nparam + bmap->n_in + bmap->n_out;
 831         total = total_var + bmap->n_div;
 832
 833         last_var = total - 1;
 834         for (done = 0; done < bmap->n_eq; ++done) {
 835                 for (; last_var >= 0; --last_var) {
 836                         for (k = done; k < bmap->n_eq; ++k)
 837                                 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
 838                                         break;
 839                         if (k < bmap->n_eq)
 840                                 break;
 841                 }
 842                 if (last_var < 0)
 843                         break;
 844                 if (k != done)
 845                         swap_equality(bmap, k, done);
 846                 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
 847                         isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
 848
 849                 for (k = 0; k < bmap->n_eq; ++k) {
 850                         if (k == done)
 851                                 continue;
 852                         if (isl_int_is_zero(bmap->eq[k][1+last_var]))
 853                                 continue;
 854                         if (progress)
 855                                 *progress = 1;
 856                         isl_seq_elim(bmap->eq[k], bmap->eq[done],
 857                                         1+last_var, 1+total, NULL);
 858                 }
 859
 860                 for (k = 0; k < bmap->n_ineq; ++k) {
 861                         if (isl_int_is_zero(bmap->ineq[k][1+last_var]))
 862                                 continue;
 863                         if (progress)
 864                                 *progress = 1;
 865                         isl_seq_elim(bmap->ineq[k], bmap->eq[done],
 866                                         1+last_var, 1+total, NULL);
 867                 }
 868
 869                 for (k = 0; k < bmap->n_div; ++k) {
 870                         if (isl_int_is_zero(bmap->div[k][0]))
 871                                 continue;
 872                         if (isl_int_is_zero(bmap->div[k][1+1+last_var]))
 873                                 continue;
 874                         if (progress)
 875                                 *progress = 1;
 876                         isl_seq_elim(bmap->div[k]+1, bmap->eq[done],
 877                                         1+last_var, 1+total, &bmap->div[k][0]);
 878                 }
 879
 880                 if (last_var >= total_var &&
 881                     isl_int_is_zero(bmap->div[last_var - total_var][0])) {
 882                         unsigned div = last_var - total_var;
 883                         isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
 884                         isl_int_set_si(bmap->div[div][1+1+last_var], 0);
 885                         isl_int_set(bmap->div[div][0],
 886                                     bmap->eq[done][1+last_var]);
 887                 }
 888         }
 889         if (done == bmap->n_eq)
 890                 return bmap;
 891         for (k = done; k < bmap->n_eq; ++k) {
 892                 if (isl_int_is_zero(bmap->eq[k][0]))
 893                         continue;
 894                 return isl_basic_map_set_to_empty(bmap);
 895         }
 896         isl_basic_map_free_equality(bmap, bmap->n_eq-done);
 897         return bmap;
 898 }
 899
 900 struct isl_basic_set *isl_basic_set_gauss(
 901         struct isl_basic_set *bset, int *progress)
 902 {
 903         return (struct isl_basic_set*)isl_basic_map_gauss(
 904                         (struct isl_basic_map *)bset, progress);
 905 }
 906
 907 static unsigned int round_up(unsigned int v)
 908 {
 909         int old_v = v;
 910
 911         while (v) {
 912                 old_v = v;
 913                 v ^= v & -v;
 914         }
 915         return old_v << 1;
 916 }
 917
 918 static int hash_index(int *index, unsigned int size, int bits,
 919                         struct isl_basic_map *bmap, int k)
 920 {
 921         int h;
 922         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 923         u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
 924         for (h = hash; index[h]; h = (h+1) % size)
 925                 if (k != index[h]-1 &&
 926                     isl_seq_eq(bmap->ineq[k]+1,
 927                                bmap->ineq[index[h]-1]+1, total))
 928                         break;
 929         return h;
 930 }
 931
 932 static struct isl_basic_map *remove_duplicate_constraints(
 933         struct isl_basic_map *bmap, int *progress)
 934 {
 935         unsigned int size;
 936         int *index;
 937         int k, l, h;
 938         int bits;
 939         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 940         isl_int sum;
 941
 942         if (bmap->n_ineq <= 1)
 943                 return bmap;
 944
 945         size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
 946         bits = ffs(size) - 1;
 947         index = isl_alloc_array(ctx, int, size);
 948         memset(index, 0, size * sizeof(int));
 949         if (!index)
 950                 return bmap;
 951
 952         index[isl_seq_hash(bmap->ineq[0]+1, total, bits)] = 1;
 953         for (k = 1; k < bmap->n_ineq; ++k) {
 954                 h = hash_index(index, size, bits, bmap, k);
 955                 if (!index[h]) {
 956                         index[h] = k+1;
 957                         continue;
 958                 }
 959                 if (progress)
 960                         *progress = 1;
 961                 l = index[h] - 1;
 962                 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
 963                         swap_inequality(bmap, k, l);
 964                 isl_basic_map_drop_inequality(bmap, k);
 965                 --k;
 966         }
 967         isl_int_init(sum);
 968         for (k = 0; k < bmap->n_ineq-1; ++k) {
 969                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
 970                 h = hash_index(index, size, bits, bmap, k);
 971                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
 972                 if (!index[h])
 973                         continue;
 974                 l = index[h] - 1;
 975                 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
 976                 if (isl_int_is_pos(sum))
 977                         continue;
 978                 if (isl_int_is_zero(sum)) {
 979                         /* We need to break out of the loop after these
 980                          * changes since the contents of the hash
 981                          * will no longer be valid.
 982                          * Plus, we probably we want to regauss first.
 983                          */
 984                         isl_basic_map_drop_inequality(bmap, l);
 985                         isl_basic_map_inequality_to_equality(bmap, k);
 986                 } else
 987                         bmap = isl_basic_map_set_to_empty(bmap);
 988                 break;
 989         }
 990         isl_int_clear(sum);
 991
 992         free(index);
 993         return bmap;
 994 }
 995
 996 static struct isl_basic_map *remove_duplicate_divs(
 997         struct isl_basic_map *bmap, int *progress)
 998 {
 999         unsigned int size;
1000         int *index;
1001         int k, l, h;
1002         int bits;
1003         struct isl_blk eq;
1004         unsigned total_var = bmap->nparam + bmap->n_in + bmap->n_out;
1005         unsigned total = total_var + bmap->n_div;
1006         struct isl_ctx *ctx;
1007
1008         if (bmap->n_div <= 1)
1009                 return bmap;
1010
1011         ctx = bmap->ctx;
1012         for (k = bmap->n_div - 1; k >= 0; --k)
1013                 if (!isl_int_is_zero(bmap->div[k][0]))
1014                         break;
1015         if (k <= 0)
1016                 return bmap;
1017
1018         size = round_up(4 * bmap->n_div / 3 - 1);
1019         bits = ffs(size) - 1;
1020         index = isl_alloc_array(ctx, int, size);
1021         memset(index, 0, size * sizeof(int));
1022         if (!index)
1023                 return bmap;
1024         eq = isl_blk_alloc(ctx, 1+total);
1025         if (isl_blk_is_error(eq))
1026                 goto out;
1027
1028         isl_seq_clr(eq.data, 1+total);
1029         index[isl_seq_hash(bmap->div[k], 2+total, bits)] = k + 1;
1030         for (--k; k >= 0; --k) {
1031                 u_int32_t hash;
1032
1033                 if (isl_int_is_zero(bmap->div[k][0]))
1034                         continue;
1035
1036                 hash = isl_seq_hash(bmap->div[k], 2+total, bits);
1037                 for (h = hash; index[h]; h = (h+1) % size)
1038                         if (isl_seq_eq(bmap->div[k],
1039                                        bmap->div[index[h]-1], 2+total))
1040                                 break;
1041                 if (index[h]) {
1042                         *progress = 1;
1043                         l = index[h] - 1;
1044                         isl_int_set_si(eq.data[1+total_var+k], -1);
1045                         isl_int_set_si(eq.data[1+total_var+l], 1);
1046                         eliminate_div(bmap, eq.data, l);
1047                         isl_int_set_si(eq.data[1+total_var+k], 0);
1048                         isl_int_set_si(eq.data[1+total_var+l], 0);
1049                 }
1050                 index[h] = k+1;
1051         }
1052
1053         isl_blk_free(ctx, eq);
1054 out:
1055         free(index);
1056         return bmap;
1057 }
1058
1059 /* Elimininate divs based on equalities
1060  */
1061 static struct isl_basic_map *eliminate_divs_eq(
1062                 struct isl_basic_map *bmap, int *progress)
1063 {
1064         int d;
1065         int i;
1066         int modified = 0;
1067         unsigned off;
1068
1069         if (!bmap)
1070                 return NULL;
1071
1072         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1073
1074         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1075                 for (i = 0; i < bmap->n_eq; ++i) {
1076                         if (!isl_int_is_one(bmap->eq[i][off + d]) &&
1077                             !isl_int_is_negone(bmap->eq[i][off + d]))
1078                                 continue;
1079                         modified = 1;
1080                         *progress = 1;
1081                         eliminate_div(bmap, bmap->eq[i], d);
1082                         isl_basic_map_drop_equality(bmap, i);
1083                         break;
1084                 }
1085         }
1086         if (modified)
1087                 return eliminate_divs_eq(bmap, progress);
1088         return bmap;
1089 }
1090
1091 static struct isl_basic_map *normalize_constraints(struct isl_basic_map *bmap)
1092 {
1093         int i;
1094         isl_int gcd;
1095         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1096
1097         isl_int_init(gcd);
1098         for (i = bmap->n_eq - 1; i >= 0; --i) {
1099                 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
1100                 if (isl_int_is_zero(gcd)) {
1101                         if (!isl_int_is_zero(bmap->eq[i][0])) {
1102                                 bmap = isl_basic_map_set_to_empty(bmap);
1103                                 break;
1104                         }
1105                         isl_basic_map_drop_equality(bmap, i);
1106                         continue;
1107                 }
1108                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1109                         isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
1110                 if (isl_int_is_one(gcd))
1111                         continue;
1112                 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
1113                         bmap = isl_basic_map_set_to_empty(bmap);
1114                         break;
1115                 }
1116                 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
1117         }
1118
1119         for (i = bmap->n_ineq - 1; i >= 0; --i) {
1120                 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
1121                 if (isl_int_is_zero(gcd)) {
1122                         if (isl_int_is_neg(bmap->ineq[i][0])) {
1123                                 bmap = isl_basic_map_set_to_empty(bmap);
1124                                 break;
1125                         }
1126                         isl_basic_map_drop_inequality(bmap, i);
1127                         continue;
1128                 }
1129                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1130                         isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
1131                 if (isl_int_is_one(gcd))
1132                         continue;
1133                 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
1134                 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
1135         }
1136         isl_int_clear(gcd);
1137
1138         return bmap;
1139 }
1140
1141
1142 /* Eliminate the specified variables from the constraints using
1143  * Fourier-Motzkin.  The variables themselves are not removed.
1144  */
1145 struct isl_basic_map *isl_basic_map_eliminate_vars(
1146         struct isl_basic_map *bmap, int pos, unsigned n)
1147 {
1148         int d;
1149         int i, j, k;
1150         unsigned total;
1151         unsigned extra;
1152
1153         if (!bmap)
1154                 return NULL;
1155         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1156         extra = bmap->extra - bmap->n_div;
1157
1158         bmap = isl_basic_map_cow(bmap);
1159         bmap = isl_basic_map_gauss(bmap, NULL);
1160         for (d = pos + n - 1; d >= pos; --d) {
1161                 int n_lower, n_upper;
1162                 if (!bmap)
1163                         return NULL;
1164                 for (i = 0; i < bmap->n_eq; ++i) {
1165                         if (isl_int_is_zero(bmap->eq[i][1+d]))
1166                                 continue;
1167                         isl_basic_map_drop_equality(bmap, i);
1168                         break;
1169                 }
1170                 if (i < bmap->n_eq)
1171                         continue;
1172                 n_lower = 0;
1173                 n_upper = 0;
1174                 for (i = 0; i < bmap->n_ineq; ++i) {
1175                         if (isl_int_is_pos(bmap->ineq[i][1+d]))
1176                                 n_lower++;
1177                         else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1178                                 n_upper++;
1179                 }
1180                 bmap = isl_basic_map_extend(bmap,
1181                                 bmap->nparam, bmap->n_in, bmap->n_out, 0,
1182                                 0, n_lower * n_upper);
1183                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1184                         int last;
1185                         if (isl_int_is_zero(bmap->ineq[i][1+d]))
1186                                 continue;
1187                         last = -1;
1188                         for (j = 0; j < i; ++j) {
1189                                 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1190                                         continue;
1191                                 last = j;
1192                                 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1193                                     isl_int_sgn(bmap->ineq[j][1+d]))
1194                                         continue;
1195                                 k = isl_basic_map_alloc_inequality(bmap);
1196                                 if (k < 0)
1197                                         goto error;
1198                                 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1199                                                 1+total);
1200                                 isl_seq_clr(bmap->ineq[k]+1+total, extra);
1201                                 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1202                                                 1+d, 1+total, NULL);
1203                         }
1204                         isl_basic_map_drop_inequality(bmap, i);
1205                         i = last + 1;
1206                 }
1207                 if (n_lower > 0 && n_upper > 0) {
1208                         bmap = normalize_constraints(bmap);
1209                         bmap = remove_duplicate_constraints(bmap, NULL);
1210                         bmap = isl_basic_map_gauss(bmap, NULL);
1211                         bmap = isl_basic_map_convex_hull(bmap);
1212                         if (!bmap)
1213                                 goto error;
1214                         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1215                                 break;
1216                 }
1217         }
1218         return bmap;
1219 error:
1220         isl_basic_map_free(bmap);
1221         return NULL;
1222 }
1223
1224 struct isl_basic_set *isl_basic_set_eliminate_vars(
1225         struct isl_basic_set *bset, unsigned pos, unsigned n)
1226 {
1227         return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1228                         (struct isl_basic_map *)bset, pos, n);
1229 }
1230
1231 /* Project out n dimensions starting at first using Fourier-Motzkin */
1232 struct isl_set *isl_set_remove_dims(struct isl_set *set,
1233         unsigned first, unsigned n)
1234 {
1235         int i;
1236
1237         if (n == 0)
1238                 return set;
1239
1240         set = isl_set_cow(set);
1241         if (!set)
1242                 return NULL;
1243         isl_assert(set->ctx, first+n <= set->dim, goto error);
1244
1245         for (i = 0; i < set->n; ++i) {
1246                 set->p[i] = isl_basic_set_eliminate_vars(set->p[i],
1247                                                         set->nparam + first, n);
1248                 if (!set->p[i])
1249                         goto error;
1250         }
1251         set = isl_set_drop_dims(set, first, n);
1252         return set;
1253 error:
1254         isl_set_free(set);
1255         return NULL;
1256 }
1257
1258 /* Project out n dimensions starting at first using Fourier-Motzkin */
1259 struct isl_basic_set *isl_basic_set_remove_dims(struct isl_basic_set *bset,
1260         unsigned first, unsigned n)
1261 {
1262         bset = isl_basic_set_eliminate_vars(bset, bset->nparam + first, n);
1263         bset = isl_basic_set_drop_dims(bset, first, n);
1264         return bset;
1265 }
1266
1267 /* Elimininate divs based on inequalities
1268  */
1269 static struct isl_basic_map *eliminate_divs_ineq(
1270                 struct isl_basic_map *bmap, int *progress)
1271 {
1272         int d;
1273         int i;
1274         unsigned off;
1275         struct isl_ctx *ctx;
1276
1277         if (!bmap)
1278                 return NULL;
1279
1280         ctx = bmap->ctx;
1281         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1282
1283         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1284                 for (i = 0; i < bmap->n_eq; ++i)
1285                         if (!isl_int_is_zero(bmap->eq[i][off + d]))
1286                                 break;
1287                 if (i < bmap->n_eq)
1288                         continue;
1289                 for (i = 0; i < bmap->n_ineq; ++i)
1290                         if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
1291                                 break;
1292                 if (i < bmap->n_ineq)
1293                         continue;
1294                 *progress = 1;
1295                 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
1296                 if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1297                         break;
1298                 bmap = isl_basic_map_drop_div(bmap, d);
1299                 if (!bmap)
1300                         break;
1301         }
1302         return bmap;
1303 }
1304
1305 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1306 {
1307         int progress = 1;
1308         if (!bmap)
1309                 return NULL;
1310         while (progress) {
1311                 progress = 0;
1312                 bmap = normalize_constraints(bmap);
1313                 bmap = eliminate_divs_eq(bmap, &progress);
1314                 bmap = eliminate_divs_ineq(bmap, &progress);
1315                 bmap = isl_basic_map_gauss(bmap, &progress);
1316                 bmap = remove_duplicate_divs(bmap, &progress);
1317                 bmap = remove_duplicate_constraints(bmap, &progress);
1318         }
1319         return bmap;
1320 }
1321
1322 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1323 {
1324         return (struct isl_basic_set *)
1325                 isl_basic_map_simplify((struct isl_basic_map *)bset);
1326 }
1327
1328 static void dump_term(struct isl_basic_map *bmap,
1329                         isl_int c, int pos, FILE *out)
1330 {
1331         unsigned in = bmap->n_in;
1332         unsigned dim = bmap->n_in + bmap->n_out;
1333         if (!pos)
1334                 isl_int_print(out, c);
1335         else {
1336                 if (!isl_int_is_one(c))
1337                         isl_int_print(out, c);
1338                 if (pos < 1 + bmap->nparam)
1339                         fprintf(out, "p%d", pos - 1);
1340                 else if (pos < 1 + bmap->nparam + in)
1341                         fprintf(out, "i%d", pos - 1 - bmap->nparam);
1342                 else if (pos < 1 + bmap->nparam + dim)
1343                         fprintf(out, "o%d", pos - 1 - bmap->nparam - in);
1344                 else
1345                         fprintf(out, "e%d", pos - 1 - bmap->nparam - dim);
1346         }
1347 }
1348
1349 static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
1350                                 int sign, FILE *out)
1351 {
1352         int i;
1353         int first;
1354         unsigned len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1355         isl_int v;
1356
1357         isl_int_init(v);
1358         for (i = 0, first = 1; i < len; ++i) {
1359                 if (isl_int_sgn(c[i]) * sign <= 0)
1360                         continue;
1361                 if (!first)
1362                         fprintf(out, " + ");
1363                 first = 0;
1364                 isl_int_abs(v, c[i]);
1365                 dump_term(bmap, v, i, out);
1366         }
1367         isl_int_clear(v);
1368         if (first)
1369                 fprintf(out, "0");
1370 }
1371
1372 static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
1373                                 const char *op, FILE *out, int indent)
1374 {
1375         int i;
1376
1377         fprintf(out, "%*s", indent, "");
1378
1379         dump_constraint_sign(bmap, c, 1, out);
1380         fprintf(out, " %s ", op);
1381         dump_constraint_sign(bmap, c, -1, out);
1382
1383         fprintf(out, "\n");
1384
1385         for (i = bmap->n_div; i < bmap->extra; ++i) {
1386                 if (isl_int_is_zero(c[1+bmap->nparam+bmap->n_in+bmap->n_out+i]))
1387                         continue;
1388                 fprintf(out, "%*s", indent, "");
1389                 fprintf(out, "ERROR: unused div coefficient not zero\n");
1390         }
1391 }
1392
1393 static void dump_constraints(struct isl_basic_map *bmap,
1394                                 isl_int **c, unsigned n,
1395                                 const char *op, FILE *out, int indent)
1396 {
1397         int i;
1398
1399         for (i = 0; i < n; ++i)
1400                 dump_constraint(bmap, c[i], op, out, indent);
1401 }
1402
1403 static void dump_affine(struct isl_basic_map *bmap, isl_int *exp, FILE *out)
1404 {
1405         int j;
1406         int first = 1;
1407         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1408
1409         for (j = 0; j < 1 + total; ++j) {
1410                 if (isl_int_is_zero(exp[j]))
1411                         continue;
1412                 if (!first && isl_int_is_pos(exp[j]))
1413                         fprintf(out, "+");
1414                 dump_term(bmap, exp[j], j, out);
1415                 first = 0;
1416         }
1417 }
1418
1419 static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
1420 {
1421         int i;
1422
1423         dump_constraints(bmap, bmap->eq, bmap->n_eq, "=", out, indent);
1424         dump_constraints(bmap, bmap->ineq, bmap->n_ineq, ">=", out, indent);
1425
1426         for (i = 0; i < bmap->n_div; ++i) {
1427                 fprintf(out, "%*s", indent, "");
1428                 fprintf(out, "e%d = [(", i);
1429                 dump_affine(bmap, bmap->div[i]+1, out);
1430                 fprintf(out, ")/");
1431                 isl_int_print(out, bmap->div[i][0]);
1432                 fprintf(out, "]\n");
1433         }
1434 }
1435
1436 void isl_basic_set_dump(struct isl_basic_set *bset, FILE *out, int indent)
1437 {
1438         if (!bset) {
1439                 fprintf(out, "null basic set\n");
1440                 return;
1441         }
1442
1443         fprintf(out, "%*s", indent, "");
1444         fprintf(out, "ref: %d, nparam: %d, dim: %d, extra: %d, flags: %x\n",
1445                         bset->ref, bset->nparam, bset->dim, bset->extra,
1446                         bset->flags);
1447         dump((struct isl_basic_map *)bset, out, indent);
1448 }
1449
1450 void isl_basic_map_dump(struct isl_basic_map *bmap, FILE *out, int indent)
1451 {
1452         if (!bmap) {
1453                 fprintf(out, "null basic map\n");
1454                 return;
1455         }
1456
1457         fprintf(out, "%*s", indent, "");
1458         fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d, flags: %x\n",
1459                 bmap->ref,
1460                 bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, bmap->flags);
1461         dump(bmap, out, indent);
1462 }
1463
1464 int isl_inequality_negate(struct isl_basic_map *bmap, unsigned pos)
1465 {
1466         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1467         if (!bmap)
1468                 return -1;
1469         isl_assert(bmap->ctx, pos < bmap->n_ineq, return -1);
1470         isl_seq_neg(bmap->ineq[pos], bmap->ineq[pos], 1 + total);
1471         isl_int_sub_ui(bmap->ineq[pos][0], bmap->ineq[pos][0], 1);
1472         return 0;
1473 }
1474
1475 struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
1476                 unsigned nparam, unsigned dim, int n, unsigned flags)
1477 {
1478         struct isl_set *set;
1479
1480         isl_assert(ctx, n >= 0, return NULL);
1481         set = isl_alloc(ctx, struct isl_set,
1482                         sizeof(struct isl_set) +
1483                         n * sizeof(struct isl_basic_set *));
1484         if (!set)
1485                 return NULL;
1486
1487         set->ctx = ctx;
1488         isl_ctx_ref(ctx);
1489         set->ref = 1;
1490         set->size = n;
1491         set->n = 0;
1492         set->nparam = nparam;
1493         set->zero = 0;
1494         set->dim = dim;
1495         set->flags = flags;
1496         return set;
1497 }
1498
1499 struct isl_set *isl_set_dup(struct isl_set *set)
1500 {
1501         int i;
1502         struct isl_set *dup;
1503
1504         if (!set)
1505                 return NULL;
1506
1507         dup = isl_set_alloc(set->ctx, set->nparam, set->dim, set->n, set->flags);
1508         if (!dup)
1509                 return NULL;
1510         for (i = 0; i < set->n; ++i)
1511                 dup = isl_set_add(dup,
1512                                   isl_basic_set_dup(set->p[i]));
1513         return dup;
1514 }
1515
1516 struct isl_set *isl_set_from_basic_set(struct isl_basic_set *bset)
1517 {
1518         struct isl_set *set;
1519
1520         if (!bset)
1521                 return NULL;
1522
1523         set = isl_set_alloc(bset->ctx, bset->nparam, bset->dim, 1, ISL_MAP_DISJOINT);
1524         if (!set) {
1525                 isl_basic_set_free(bset);
1526                 return NULL;
1527         }
1528         return isl_set_add(set, bset);
1529 }
1530
1531 struct isl_map *isl_map_from_basic_map(struct isl_basic_map *bmap)
1532 {
1533         struct isl_map *map;
1534
1535         if (!bmap)
1536                 return NULL;
1537
1538         map = isl_map_alloc(bmap->ctx, bmap->nparam, bmap->n_in, bmap->n_out, 1,
1539                                 ISL_MAP_DISJOINT);
1540         if (!map) {
1541                 isl_basic_map_free(bmap);
1542                 return NULL;
1543         }
1544         return isl_map_add(map, bmap);
1545 }
1546
1547 struct isl_set *isl_set_add(struct isl_set *set, struct isl_basic_set *bset)
1548 {
1549         if (!bset || !set)
1550                 goto error;
1551         isl_assert(set->ctx, set->nparam == bset->nparam, goto error);
1552         isl_assert(set->ctx, set->dim == bset->dim, goto error);
1553         isl_assert(set->ctx, set->n < set->size, goto error);
1554         set->p[set->n] = bset;
1555         set->n++;
1556         return set;
1557 error:
1558         if (set)
1559                 isl_set_free(set);
1560         if (bset)
1561                 isl_basic_set_free(bset);
1562         return NULL;
1563 }
1564
1565 void isl_set_free(struct isl_set *set)
1566 {
1567         int i;
1568
1569         if (!set)
1570                 return;
1571
1572         if (--set->ref > 0)
1573                 return;
1574
1575         isl_ctx_deref(set->ctx);
1576         for (i = 0; i < set->n; ++i)
1577                 isl_basic_set_free(set->p[i]);
1578         free(set);
1579 }
1580
1581 void isl_set_dump(struct isl_set *set, FILE *out, int indent)
1582 {
1583         int i;
1584
1585         if (!set) {
1586                 fprintf(out, "null set\n");
1587                 return;
1588         }
1589
1590         fprintf(out, "%*s", indent, "");
1591         fprintf(out, "ref: %d, n: %d, nparam: %d, dim: %d, flags: %x\n",
1592                         set->ref, set->n, set->nparam, set->dim, set->flags);
1593         for (i = 0; i < set->n; ++i) {
1594                 fprintf(out, "%*s", indent, "");
1595                 fprintf(out, "basic set %d:\n", i);
1596                 isl_basic_set_dump(set->p[i], out, indent+4);
1597         }
1598 }
1599
1600 void isl_map_dump(struct isl_map *map, FILE *out, int indent)
1601 {
1602         int i;
1603
1604         if (!map) {
1605                 fprintf(out, "null map\n");
1606                 return;
1607         }
1608
1609         fprintf(out, "%*s", indent, "");
1610         fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d, flags: %x\n",
1611                         map->ref, map->n, map->nparam, map->n_in, map->n_out,
1612                         map->flags);
1613         for (i = 0; i < map->n; ++i) {
1614                 fprintf(out, "%*s", indent, "");
1615                 fprintf(out, "basic map %d:\n", i);
1616                 isl_basic_map_dump(map->p[i], out, indent+4);
1617         }
1618 }
1619
1620 struct isl_basic_map *isl_basic_map_intersect_domain(
1621                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1622 {
1623         struct isl_basic_map *bmap_domain;
1624
1625         if (!bmap || !bset)
1626                 goto error;
1627
1628         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
1629         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1630
1631         bmap = isl_basic_map_extend(bmap,
1632                         bset->nparam, bmap->n_in, bmap->n_out,
1633                         bset->n_div, bset->n_eq, bset->n_ineq);
1634         if (!bmap)
1635                 goto error;
1636         bmap_domain = isl_basic_map_from_basic_set(bset,
1637                                                 bset->dim, 0);
1638         bmap = add_constraints(bmap, bmap_domain, 0, 0);
1639
1640         bmap = isl_basic_map_simplify(bmap);
1641         return isl_basic_map_finalize(bmap);
1642 error:
1643         isl_basic_map_free(bmap);
1644         isl_basic_set_free(bset);
1645         return NULL;
1646 }
1647
1648 struct isl_basic_map *isl_basic_map_intersect_range(
1649                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1650 {
1651         struct isl_basic_map *bmap_range;
1652
1653         if (!bmap || !bset)
1654                 goto error;
1655
1656         isl_assert(bset->ctx, bset->dim == bmap->n_out, goto error);
1657         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1658
1659         bmap = isl_basic_map_extend(bmap,
1660                         bset->nparam, bmap->n_in, bmap->n_out,
1661                         bset->n_div, bset->n_eq, bset->n_ineq);
1662         if (!bmap)
1663                 goto error;
1664         bmap_range = isl_basic_map_from_basic_set(bset,
1665                                                 0, bset->dim);
1666         bmap = add_constraints(bmap, bmap_range, 0, 0);
1667
1668         bmap = isl_basic_map_simplify(bmap);
1669         return isl_basic_map_finalize(bmap);
1670 error:
1671         isl_basic_map_free(bmap);
1672         isl_basic_set_free(bset);
1673         return NULL;
1674 }
1675
1676 struct isl_basic_map *isl_basic_map_intersect(
1677                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1678 {
1679         if (!bmap1 || !bmap2)
1680                 goto error;
1681
1682         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
1683         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
1684         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
1685
1686         bmap1 = isl_basic_map_extend(bmap1,
1687                         bmap1->nparam, bmap1->n_in, bmap1->n_out,
1688                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
1689         if (!bmap1)
1690                 goto error;
1691         bmap1 = add_constraints(bmap1, bmap2, 0, 0);
1692
1693         bmap1 = isl_basic_map_simplify(bmap1);
1694         return isl_basic_map_finalize(bmap1);
1695 error:
1696         isl_basic_map_free(bmap1);
1697         isl_basic_map_free(bmap2);
1698         return NULL;
1699 }
1700
1701 struct isl_basic_set *isl_basic_set_intersect(
1702                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1703 {
1704         return (struct isl_basic_set *)
1705                 isl_basic_map_intersect(
1706                         (struct isl_basic_map *)bset1,
1707                         (struct isl_basic_map *)bset2);
1708 }
1709
1710 struct isl_map *isl_map_intersect(struct isl_map *map1, struct isl_map *map2)
1711 {
1712         unsigned flags = 0;
1713         struct isl_map *result;
1714         int i, j;
1715
1716         if (!map1 || !map2)
1717                 goto error;
1718
1719         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
1720             F_ISSET(map2, ISL_MAP_DISJOINT))
1721                 FL_SET(flags, ISL_MAP_DISJOINT);
1722
1723         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
1724                                 map1->n * map2->n, flags);
1725         if (!result)
1726                 goto error;
1727         for (i = 0; i < map1->n; ++i)
1728                 for (j = 0; j < map2->n; ++j) {
1729                         struct isl_basic_map *part;
1730                         part = isl_basic_map_intersect(
1731                                     isl_basic_map_copy(map1->p[i]),
1732                                     isl_basic_map_copy(map2->p[j]));
1733                         if (isl_basic_map_is_empty(part))
1734                                 isl_basic_map_free(part);
1735                         else
1736                                 result = isl_map_add(result, part);
1737                         if (!result)
1738                                 goto error;
1739                 }
1740         isl_map_free(map1);
1741         isl_map_free(map2);
1742         return result;
1743 error:
1744         isl_map_free(map1);
1745         isl_map_free(map2);
1746         return NULL;
1747 }
1748
1749 struct isl_set *isl_set_intersect(struct isl_set *set1, struct isl_set *set2)
1750 {
1751         return (struct isl_set *)
1752                 isl_map_intersect((struct isl_map *)set1,
1753                                   (struct isl_map *)set2);
1754 }
1755
1756 struct isl_basic_map *isl_basic_map_reverse(struct isl_basic_map *bmap)
1757 {
1758         struct isl_basic_set *bset;
1759         unsigned in;
1760
1761         if (!bmap)
1762                 return NULL;
1763         bmap = isl_basic_map_cow(bmap);
1764         if (!bmap)
1765                 return NULL;
1766         in = bmap->n_in;
1767         bset = isl_basic_set_from_basic_map(bmap);
1768         bset = isl_basic_set_swap_vars(bset, in);
1769         return isl_basic_map_from_basic_set(bset, bset->dim-in, in);
1770 }
1771
1772 /* Turn final n dimensions into existentially quantified variables.
1773  */
1774 struct isl_basic_set *isl_basic_set_project_out(
1775                 struct isl_basic_set *bset, unsigned n, unsigned flags)
1776 {
1777         int i;
1778         size_t row_size;
1779         isl_int **new_div;
1780         isl_int *old;
1781
1782         if (!bset)
1783                 return NULL;
1784
1785         isl_assert(bset->ctx, n <= bset->dim, goto error);
1786
1787         if (n == 0)
1788                 return bset;
1789
1790         bset = isl_basic_set_cow(bset);
1791
1792         row_size = 1 + bset->nparam + bset->dim + bset->extra;
1793         old = bset->block2.data;
1794         bset->block2 = isl_blk_extend(bset->ctx, bset->block2,
1795                                         (bset->extra + n) * (1 + row_size));
1796         if (!bset->block2.data)
1797                 goto error;
1798         new_div = isl_alloc_array(ctx, isl_int *, bset->extra + n);
1799         if (!new_div)
1800                 goto error;
1801         for (i = 0; i < n; ++i) {
1802                 new_div[i] = bset->block2.data +
1803                                 (bset->extra + i) * (1 + row_size);
1804                 isl_seq_clr(new_div[i], 1 + row_size);
1805         }
1806         for (i = 0; i < bset->extra; ++i)
1807                 new_div[n + i] = bset->block2.data + (bset->div[i] - old);
1808         free(bset->div);
1809         bset->div = new_div;
1810         bset->n_div += n;
1811         bset->extra += n;
1812         bset->dim -= n;
1813         bset = isl_basic_set_simplify(bset);
1814         return isl_basic_set_finalize(bset);
1815 error:
1816         isl_basic_set_free(bset);
1817         return NULL;
1818 }
1819
1820 struct isl_basic_map *isl_basic_map_apply_range(
1821                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1822 {
1823         struct isl_basic_set *bset;
1824         unsigned n_in, n_out;
1825
1826         if (!bmap1 || !bmap2)
1827                 goto error;
1828
1829         isl_assert(bmap->ctx, bmap1->n_out == bmap2->n_in, goto error);
1830         isl_assert(bmap->ctx, bmap1->nparam == bmap2->nparam, goto error);
1831
1832         n_in = bmap1->n_in;
1833         n_out = bmap2->n_out;
1834
1835         bmap2 = isl_basic_map_reverse(bmap2);
1836         if (!bmap2)
1837                 goto error;
1838         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
1839                         bmap1->n_in + bmap2->n_in, bmap1->n_out,
1840                         bmap2->extra, bmap2->n_eq, bmap2->n_ineq);
1841         if (!bmap1)
1842                 goto error;
1843         bmap1 = add_constraints(bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
1844         bmap1 = isl_basic_map_simplify(bmap1);
1845         bset = isl_basic_set_from_basic_map(bmap1);
1846         bset = isl_basic_set_project_out(bset,
1847                                                 bset->dim - (n_in + n_out), 0);
1848         return isl_basic_map_from_basic_set(bset, n_in, n_out);
1849 error:
1850         isl_basic_map_free(bmap1);
1851         isl_basic_map_free(bmap2);
1852         return NULL;
1853 }
1854
1855 struct isl_basic_set *isl_basic_set_apply(
1856                 struct isl_basic_set *bset, struct isl_basic_map *bmap)
1857 {
1858         if (!bset || !bmap)
1859                 goto error;
1860
1861         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
1862
1863         return (struct isl_basic_set *)
1864                 isl_basic_map_apply_range((struct isl_basic_map *)bset, bmap);
1865 error:
1866         isl_basic_set_free(bset);
1867         isl_basic_map_free(bmap);
1868         return NULL;
1869 }
1870
1871 struct isl_basic_map *isl_basic_map_apply_domain(
1872                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1873 {
1874         if (!bmap1 || !bmap2)
1875                 goto error;
1876
1877         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
1878         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
1879
1880         bmap1 = isl_basic_map_reverse(bmap1);
1881         bmap1 = isl_basic_map_apply_range(bmap1, bmap2);
1882         return isl_basic_map_reverse(bmap1);
1883 error:
1884         isl_basic_map_free(bmap1);
1885         isl_basic_map_free(bmap2);
1886         return NULL;
1887 }
1888
1889 static struct isl_basic_map *var_equal(struct isl_ctx *ctx,
1890                 struct isl_basic_map *bmap, unsigned pos)
1891 {
1892         int i;
1893         i = isl_basic_map_alloc_equality(bmap);
1894         if (i < 0)
1895                 goto error;
1896         isl_seq_clr(bmap->eq[i],
1897                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1898         isl_int_set_si(bmap->eq[i][1+bmap->nparam+pos], -1);
1899         isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1900         return isl_basic_map_finalize(bmap);
1901 error:
1902         isl_basic_map_free(bmap);
1903         return NULL;
1904 }
1905
1906 static struct isl_basic_map *var_more(struct isl_ctx *ctx,
1907                 struct isl_basic_map *bmap, unsigned pos)
1908 {
1909         int i;
1910         i = isl_basic_map_alloc_inequality(bmap);
1911         if (i < 0)
1912                 goto error;
1913         isl_seq_clr(bmap->ineq[i],
1914                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1915         isl_int_set_si(bmap->ineq[i][0], -1);
1916         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], -1);
1917         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1918         return isl_basic_map_finalize(bmap);
1919 error:
1920         isl_basic_map_free(bmap);
1921         return NULL;
1922 }
1923
1924 static struct isl_basic_map *var_less(struct isl_ctx *ctx,
1925                 struct isl_basic_map *bmap, unsigned pos)
1926 {
1927         int i;
1928         i = isl_basic_map_alloc_inequality(bmap);
1929         if (i < 0)
1930                 goto error;
1931         isl_seq_clr(bmap->ineq[i],
1932                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1933         isl_int_set_si(bmap->ineq[i][0], -1);
1934         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], 1);
1935         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], -1);
1936         return isl_basic_map_finalize(bmap);
1937 error:
1938         isl_basic_map_free(bmap);
1939         return NULL;
1940 }
1941
1942 struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
1943                 unsigned nparam, unsigned in, unsigned out, unsigned n_equal)
1944 {
1945         int i;
1946         struct isl_basic_map *bmap;
1947         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, n_equal, 0);
1948         if (!bmap)
1949                 return NULL;
1950         for (i = 0; i < n_equal && bmap; ++i)
1951                 bmap = var_equal(ctx, bmap, i);
1952         return isl_basic_map_finalize(bmap);
1953 }
1954
1955 struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
1956                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1957 {
1958         int i;
1959         struct isl_basic_map *bmap;
1960         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1961         if (!bmap)
1962                 return NULL;
1963         for (i = 0; i < pos && bmap; ++i)
1964                 bmap = var_equal(ctx, bmap, i);
1965         if (bmap)
1966                 bmap = var_less(ctx, bmap, pos);
1967         return isl_basic_map_finalize(bmap);
1968 }
1969
1970 struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
1971                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1972 {
1973         int i;
1974         struct isl_basic_map *bmap;
1975         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1976         if (!bmap)
1977                 return NULL;
1978         for (i = 0; i < pos && bmap; ++i)
1979                 bmap = var_equal(ctx, bmap, i);
1980         if (bmap)
1981                 bmap = var_more(ctx, bmap, pos);
1982         return isl_basic_map_finalize(bmap);
1983 }
1984
1985 struct isl_basic_map *isl_basic_map_from_basic_set(
1986                 struct isl_basic_set *bset, unsigned n_in, unsigned n_out)
1987 {
1988         struct isl_basic_map *bmap;
1989
1990         bset = isl_basic_set_cow(bset);
1991         if (!bset)
1992                 return NULL;
1993
1994         isl_assert(ctx, bset->dim == n_in + n_out, goto error);
1995         bmap = (struct isl_basic_map *) bset;
1996         bmap->n_in = n_in;
1997         bmap->n_out = n_out;
1998         return isl_basic_map_finalize(bmap);
1999 error:
2000         isl_basic_set_free(bset);
2001         return NULL;
2002 }
2003
2004 struct isl_basic_set *isl_basic_set_from_basic_map(struct isl_basic_map *bmap)
2005 {
2006         if (!bmap)
2007                 goto error;
2008         if (bmap->n_in == 0)
2009                 return (struct isl_basic_set *)bmap;
2010         bmap = isl_basic_map_cow(bmap);
2011         if (!bmap)
2012                 goto error;
2013         bmap->n_out += bmap->n_in;
2014         bmap->n_in = 0;
2015         bmap = isl_basic_map_finalize(bmap);
2016         return (struct isl_basic_set *)bmap;
2017 error:
2018         return NULL;
2019 }
2020
2021 /* For a div d = floor(f/m), add the constraints
2022  *
2023  *              f - m d >= 0
2024  *              -(f-(n-1)) + m d >= 0
2025  *
2026  * Note that the second constraint is the negation of
2027  *
2028  *              f - m d >= n
2029  */
2030 static int add_div_constraints(struct isl_basic_map *bmap, unsigned div)
2031 {
2032         int i, j;
2033         unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2034         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
2035
2036         i = isl_basic_map_alloc_inequality(bmap);
2037         if (i < 0)
2038                 return -1;
2039         isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
2040         isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
2041
2042         j = isl_basic_map_alloc_inequality(bmap);
2043         if (j < 0)
2044                 return -1;
2045         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
2046         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
2047         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
2048         return j;
2049 }
2050
2051 struct isl_basic_set *isl_basic_map_underlying_set(
2052                 struct isl_basic_map *bmap)
2053 {
2054         if (!bmap)
2055                 goto error;
2056         if (bmap->nparam == 0 && bmap->n_in == 0 && bmap->n_div == 0)
2057                 return (struct isl_basic_set *)bmap;
2058         bmap = isl_basic_map_cow(bmap);
2059         if (!bmap)
2060                 goto error;
2061         bmap->n_out += bmap->nparam + bmap->n_in + bmap->n_div;
2062         bmap->nparam = 0;
2063         bmap->n_in = 0;
2064         bmap->extra -= bmap->n_div;
2065         bmap->n_div = 0;
2066         bmap = isl_basic_map_finalize(bmap);
2067         return (struct isl_basic_set *)bmap;
2068 error:
2069         return NULL;
2070 }
2071
2072 struct isl_basic_map *isl_basic_map_overlying_set(
2073         struct isl_basic_set *bset, struct isl_basic_map *like)
2074 {
2075         struct isl_basic_map *bmap;
2076         struct isl_ctx *ctx;
2077         unsigned total;
2078         int i, k;
2079
2080         if (!bset || !like)
2081                 goto error;
2082         ctx = bset->ctx;
2083         isl_assert(ctx, bset->dim ==
2084                         like->nparam + like->n_in + like->n_out + like->n_div,
2085                         goto error);
2086         if (like->nparam == 0 && like->n_in == 0 && like->n_div == 0) {
2087                 isl_basic_map_free(like);
2088                 return (struct isl_basic_map *)bset;
2089         }
2090         bset = isl_basic_set_cow(bset);
2091         if (!bset)
2092                 goto error;
2093         total = bset->dim + bset->extra;
2094         bmap = (struct isl_basic_map *)bset;
2095         bmap->nparam = like->nparam;
2096         bmap->n_in = like->n_in;
2097         bmap->n_out = like->n_out;
2098         bmap->extra += like->n_div;
2099         if (bmap->extra) {
2100                 unsigned ltotal;
2101                 ltotal = total - bmap->extra + like->extra;
2102                 if (ltotal > total)
2103                         ltotal = total;
2104                 bmap->block2 = isl_blk_extend(ctx, bmap->block2,
2105                                         bmap->extra * (1 + 1 + total));
2106                 if (isl_blk_is_error(bmap->block2))
2107                         goto error;
2108                 bmap->div = isl_realloc_array(ctx, bmap->div, isl_int *,
2109                                                 bmap->extra);
2110                 if (!bmap->div)
2111                         goto error;
2112                 bmap = isl_basic_map_extend(bmap, bmap->nparam,
2113                         bmap->n_in, bmap->n_out, 0, 0, 2 * like->n_div);
2114                 for (i = 0; i < like->n_div; ++i) {
2115                         k = isl_basic_map_alloc_div(bmap);
2116                         if (k < 0)
2117                                 goto error;
2118                         isl_seq_cpy(bmap->div[k], like->div[i], 1 + 1 + ltotal);
2119                         isl_seq_clr(bmap->div[k]+1+1+ltotal, total - ltotal);
2120                         if (add_div_constraints(bmap, k) < 0)
2121                                 goto error;
2122                 }
2123         }
2124         isl_basic_map_free(like);
2125         bmap = isl_basic_map_finalize(bmap);
2126         return bmap;
2127 error:
2128         isl_basic_map_free(like);
2129         isl_basic_set_free(bset);
2130         return NULL;
2131 }
2132
2133 struct isl_basic_set *isl_basic_set_from_underlying_set(
2134         struct isl_basic_set *bset, struct isl_basic_set *like)
2135 {
2136         return (struct isl_basic_set *)
2137                 isl_basic_map_overlying_set(bset, (struct isl_basic_map *)like);
2138 }
2139
2140 struct isl_set *isl_set_from_underlying_set(
2141         struct isl_set *set, struct isl_basic_set *like)
2142 {
2143         int i;
2144
2145         if (!set || !like)
2146                 goto error;
2147         isl_assert(set->ctx, set->dim == like->nparam + like->dim + like->n_div,
2148                     goto error);
2149         if (like->nparam == 0 && like->n_div == 0) {
2150                 isl_basic_set_free(like);
2151                 return set;
2152         }
2153         set = isl_set_cow(set);
2154         if (!set)
2155                 goto error;
2156         for (i = 0; i < set->n; ++i) {
2157                 set->p[i] = isl_basic_set_from_underlying_set(set->p[i],
2158                                                       isl_basic_set_copy(like));
2159                 if (!set->p[i])
2160                         goto error;
2161         }
2162         set->nparam = like->nparam;
2163         set->dim = like->dim;
2164         isl_basic_set_free(like);
2165         return set;
2166 error:
2167         isl_basic_set_free(like);
2168         isl_set_free(set);
2169         return NULL;
2170 }
2171
2172 struct isl_set *isl_map_underlying_set(struct isl_map *map)
2173 {
2174         int i;
2175
2176         map = isl_map_align_divs(map);
2177         map = isl_map_cow(map);
2178         if (!map)
2179                 return NULL;
2180
2181         for (i = 0; i < map->n; ++i) {
2182                 map->p[i] = (struct isl_basic_map *)
2183                                 isl_basic_map_underlying_set(map->p[i]);
2184                 if (!map->p[i])
2185                         goto error;
2186         }
2187         if (map->n == 0)
2188                 map->n_out += map->nparam + map->n_in;
2189         else
2190                 map->n_out = map->p[0]->n_out;
2191         map->nparam = 0;
2192         map->n_in = 0;
2193         return (struct isl_set *)map;
2194 error:
2195         isl_map_free(map);
2196         return NULL;
2197 }
2198
2199 struct isl_set *isl_set_to_underlying_set(struct isl_set *set)
2200 {
2201         return (struct isl_set *)isl_map_underlying_set((struct isl_map *)set);
2202 }
2203
2204 struct isl_basic_set *isl_basic_map_domain(struct isl_basic_map *bmap)
2205 {
2206         struct isl_basic_set *domain;
2207         unsigned n_out;
2208         if (!bmap)
2209                 return NULL;
2210         n_out = bmap->n_out;
2211         domain = isl_basic_set_from_basic_map(bmap);
2212         return isl_basic_set_project_out(domain, n_out, 0);
2213 }
2214
2215 struct isl_basic_set *isl_basic_map_range(struct isl_basic_map *bmap)
2216 {
2217         return isl_basic_map_domain(isl_basic_map_reverse(bmap));
2218 }
2219
2220 struct isl_set *isl_map_range(struct isl_map *map)
2221 {
2222         int i;
2223         struct isl_set *set;
2224
2225         if (!map)
2226                 goto error;
2227         map = isl_map_cow(map);
2228         if (!map)
2229                 goto error;
2230
2231         set = (struct isl_set *) map;
2232         set->zero = 0;
2233         for (i = 0; i < map->n; ++i) {
2234                 set->p[i] = isl_basic_map_range(map->p[i]);
2235                 if (!set->p[i])
2236                         goto error;
2237         }
2238         F_CLR(set, ISL_MAP_DISJOINT);
2239         return set;
2240 error:
2241         isl_map_free(map);
2242         return NULL;
2243 }
2244
2245 struct isl_map *isl_map_from_set(struct isl_set *set,
2246                 unsigned n_in, unsigned n_out)
2247 {
2248         int i;
2249         struct isl_map *map = NULL;
2250
2251         if (!set)
2252                 return NULL;
2253         isl_assert(set->ctx, set->dim == n_in + n_out, goto error);
2254         set = isl_set_cow(set);
2255         if (!set)
2256                 return NULL;
2257         map = (struct isl_map *)set;
2258         for (i = 0; i < set->n; ++i) {
2259                 map->p[i] = isl_basic_map_from_basic_set(
2260                                 set->p[i], n_in, n_out);
2261                 if (!map->p[i])
2262                         goto error;
2263         }
2264         map->n_in = n_in;
2265         map->n_out = n_out;
2266         return map;
2267 error:
2268         isl_set_free(set);
2269         return NULL;
2270 }
2271
2272 struct isl_set *isl_set_from_map(struct isl_map *map)
2273 {
2274         int i;
2275         struct isl_set *set = NULL;
2276
2277         if (!map)
2278                 return NULL;
2279         map = isl_map_cow(map);
2280         if (!map)
2281                 return NULL;
2282         map->n_out += map->n_in;
2283         map->n_in = 0;
2284         set = (struct isl_set *)map;
2285         for (i = 0; i < map->n; ++i) {
2286                 set->p[i] = isl_basic_set_from_basic_map(map->p[i]);
2287                 if (!set->p[i])
2288                         goto error;
2289         }
2290         return set;
2291 error:
2292         isl_map_free(map);
2293         return NULL;
2294 }
2295
2296 struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
2297                 unsigned nparam, unsigned in, unsigned out, int n,
2298                 unsigned flags)
2299 {
2300         struct isl_map *map;
2301
2302         map = isl_alloc(ctx, struct isl_map,
2303                         sizeof(struct isl_map) +
2304                         n * sizeof(struct isl_basic_map *));
2305         if (!map)
2306                 return NULL;
2307
2308         map->ctx = ctx;
2309         isl_ctx_ref(ctx);
2310         map->ref = 1;
2311         map->size = n;
2312         map->n = 0;
2313         map->nparam = nparam;
2314         map->n_in = in;
2315         map->n_out = out;
2316         map->flags = flags;
2317         return map;
2318 }
2319
2320 struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
2321                 unsigned nparam, unsigned in, unsigned out)
2322 {
2323         struct isl_basic_map *bmap;
2324         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 1, 0);
2325         bmap = isl_basic_map_set_to_empty(bmap);
2326         return bmap;
2327 }
2328
2329 struct isl_basic_set *isl_basic_set_empty(struct isl_ctx *ctx,
2330                 unsigned nparam, unsigned dim)
2331 {
2332         struct isl_basic_set *bset;
2333         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 1, 0);
2334         bset = isl_basic_set_set_to_empty(bset);
2335         return bset;
2336 }
2337
2338 struct isl_basic_map *isl_basic_map_universe(struct isl_ctx *ctx,
2339                 unsigned nparam, unsigned in, unsigned out)
2340 {
2341         struct isl_basic_map *bmap;
2342         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 0, 0);
2343         return bmap;
2344 }
2345
2346 struct isl_basic_set *isl_basic_set_universe(struct isl_ctx *ctx,
2347                 unsigned nparam, unsigned dim)
2348 {
2349         struct isl_basic_set *bset;
2350         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 0, 0);
2351         return bset;
2352 }
2353
2354 struct isl_map *isl_map_empty(struct isl_ctx *ctx,
2355                 unsigned nparam, unsigned in, unsigned out)
2356 {
2357         return isl_map_alloc(ctx, nparam, in, out, 0, ISL_MAP_DISJOINT);
2358 }
2359
2360 struct isl_set *isl_set_empty(struct isl_ctx *ctx,
2361                 unsigned nparam, unsigned dim)
2362 {
2363         return isl_set_alloc(ctx, nparam, dim, 0, ISL_MAP_DISJOINT);
2364 }
2365
2366 struct isl_map *isl_map_dup(struct isl_map *map)
2367 {
2368         int i;
2369         struct isl_map *dup;
2370
2371         if (!map)
2372                 return NULL;
2373         dup = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out, map->n,
2374                                 map->flags);
2375         for (i = 0; i < map->n; ++i)
2376                 dup = isl_map_add(dup,
2377                                   isl_basic_map_dup(map->p[i]));
2378         return dup;
2379 }
2380
2381 struct isl_map *isl_map_add(struct isl_map *map, struct isl_basic_map *bmap)
2382 {
2383         if (!bmap || !map)
2384                 goto error;
2385         isl_assert(map->ctx, map->nparam == bmap->nparam, goto error);
2386         isl_assert(map->ctx, map->n_in == bmap->n_in, goto error);
2387         isl_assert(map->ctx, map->n_out == bmap->n_out, goto error);
2388         isl_assert(map->ctx, map->n < map->size, goto error);
2389         map->p[map->n] = bmap;
2390         map->n++;
2391         return map;
2392 error:
2393         if (map)
2394                 isl_map_free(map);
2395         if (bmap)
2396                 isl_basic_map_free(bmap);
2397         return NULL;
2398 }
2399
2400 void isl_map_free(struct isl_map *map)
2401 {
2402         int i;
2403
2404         if (!map)
2405                 return;
2406
2407         if (--map->ref > 0)
2408                 return;
2409
2410         isl_ctx_deref(map->ctx);
2411         for (i = 0; i < map->n; ++i)
2412                 isl_basic_map_free(map->p[i]);
2413         free(map);
2414 }
2415
2416 struct isl_map *isl_map_extend(struct isl_map *base,
2417                 unsigned nparam, unsigned n_in, unsigned n_out)
2418 {
2419         int i;
2420
2421         base = isl_map_cow(base);
2422         if (!base)
2423                 return NULL;
2424
2425         isl_assert(base->ctx, base->nparam <= nparam, goto error);
2426         isl_assert(base->ctx, base->n_in <= n_in, goto error);
2427         isl_assert(base->ctx, base->n_out <= n_out, goto error);
2428         base->nparam = nparam;
2429         base->n_in = n_in;
2430         base->n_out = n_out;
2431         for (i = 0; i < base->n; ++i) {
2432                 base->p[i] = isl_basic_map_extend(base->p[i],
2433                                 nparam, n_in, n_out, 0, 0, 0);
2434                 if (!base->p[i])
2435                         goto error;
2436         }
2437         return base;
2438 error:
2439         isl_map_free(base);
2440         return NULL;
2441 }
2442
2443 struct isl_set *isl_set_extend(struct isl_set *base,
2444                 unsigned nparam, unsigned dim)
2445 {
2446         return (struct isl_set *)isl_map_extend((struct isl_map *)base,
2447                                                         nparam, 0, dim);
2448 }
2449
2450 static struct isl_basic_map *isl_basic_map_fix_var(struct isl_basic_map *bmap,
2451                 unsigned var, int value)
2452 {
2453         int j;
2454
2455         bmap = isl_basic_map_cow(bmap);
2456         bmap = isl_basic_map_extend(bmap,
2457                         bmap->nparam, bmap->n_in, bmap->n_out, 0, 1, 0);
2458         j = isl_basic_map_alloc_equality(bmap);
2459         if (j < 0)
2460                 goto error;
2461         isl_seq_clr(bmap->eq[j],
2462                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2463         isl_int_set_si(bmap->eq[j][1 + var], -1);
2464         isl_int_set_si(bmap->eq[j][0], value);
2465         bmap = isl_basic_map_simplify(bmap);
2466         return isl_basic_map_finalize(bmap);
2467 error:
2468         isl_basic_map_free(bmap);
2469         return NULL;
2470 }
2471
2472 struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_basic_map *bmap,
2473                 unsigned input, int value)
2474 {
2475         if (!bmap)
2476                 return NULL;
2477         isl_assert(bmap->ctx, input < bmap->n_in, goto error);
2478         return isl_basic_map_fix_var(bmap, bmap->nparam + input, value);
2479 error:
2480         isl_basic_map_free(bmap);
2481         return NULL;
2482 }
2483
2484 struct isl_basic_set *isl_basic_set_fix_dim_si(struct isl_basic_set *bset,
2485                 unsigned dim, int value)
2486 {
2487         if (!bset)
2488                 return NULL;
2489         isl_assert(bset->ctx, dim < bset->dim, goto error);
2490         return (struct isl_basic_set *)
2491                 isl_basic_map_fix_var((struct isl_basic_map *)bset,
2492                                                 bset->nparam + dim, value);
2493 error:
2494         isl_basic_set_free(bset);
2495         return NULL;
2496 }
2497
2498 struct isl_map *isl_map_fix_input_si(struct isl_map *map,
2499                 unsigned input, int value)
2500 {
2501         int i;
2502
2503         map = isl_map_cow(map);
2504         if (!map)
2505                 return NULL;
2506
2507         isl_assert(ctx, input < map->n_in, goto error);
2508         for (i = 0; i < map->n; ++i) {
2509                 map->p[i] = isl_basic_map_fix_input_si(map->p[i], input, value);
2510                 if (!map->p[i])
2511                         goto error;
2512         }
2513         return map;
2514 error:
2515         isl_map_free(map);
2516         return NULL;
2517 }
2518
2519 struct isl_set *isl_set_fix_dim_si(struct isl_set *set, unsigned dim, int value)
2520 {
2521         int i;
2522
2523         set = isl_set_cow(set);
2524         if (!set)
2525                 return NULL;
2526
2527         isl_assert(set->ctx, dim < set->dim, goto error);
2528         for (i = 0; i < set->n; ++i) {
2529                 set->p[i] = isl_basic_set_fix_dim_si(set->p[i], dim, value);
2530                 if (!set->p[i])
2531                         goto error;
2532         }
2533         return set;
2534 error:
2535         isl_set_free(set);
2536         return NULL;
2537 }
2538
2539 struct isl_map *isl_map_reverse(struct isl_map *map)
2540 {
2541         int i;
2542         unsigned t;
2543
2544         map = isl_map_cow(map);
2545         if (!map)
2546                 return NULL;
2547
2548         t = map->n_in;
2549         map->n_in = map->n_out;
2550         map->n_out = t;
2551         for (i = 0; i < map->n; ++i) {
2552                 map->p[i] = isl_basic_map_reverse(map->p[i]);
2553                 if (!map->p[i])
2554                         goto error;
2555         }
2556         return map;
2557 error:
2558         isl_map_free(map);
2559         return NULL;
2560 }
2561
2562 struct isl_map *isl_basic_map_lexmax(
2563                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2564                 struct isl_set **empty)
2565 {
2566         return isl_pip_basic_map_lexmax(bmap, dom, empty);
2567 }
2568
2569 struct isl_map *isl_basic_map_lexmin(
2570                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2571                 struct isl_set **empty)
2572 {
2573         return isl_pip_basic_map_lexmin(bmap, dom, empty);
2574 }
2575
2576 struct isl_set *isl_basic_set_lexmin(struct isl_basic_set *bset)
2577 {
2578         struct isl_basic_map *bmap = NULL;
2579         struct isl_basic_set *dom = NULL;
2580         struct isl_map *min;
2581
2582         if (!bset)
2583                 goto error;
2584         bmap = isl_basic_map_from_basic_set(bset, 0, bset->dim);
2585         if (!bmap)
2586                 goto error;
2587         dom = isl_basic_set_alloc(bmap->ctx, bmap->nparam, 0, 0, 0, 0);
2588         if (!dom)
2589                 goto error;
2590         min = isl_basic_map_lexmin(bmap, dom, NULL);
2591         return isl_map_range(min);
2592 error:
2593         isl_basic_map_free(bmap);
2594         return NULL;
2595 }
2596
2597 struct isl_map *isl_basic_map_compute_divs(struct isl_basic_map *bmap)
2598 {
2599         if (!bmap)
2600                 return NULL;
2601         if (bmap->n_div == 0)
2602                 return isl_map_from_basic_map(bmap);
2603         return isl_pip_basic_map_compute_divs(bmap);
2604 }
2605
2606 struct isl_map *isl_map_compute_divs(struct isl_map *map)
2607 {
2608         int i;
2609         struct isl_map *res;
2610
2611         if (!map)
2612                 return NULL;
2613         if (map->n == 0)
2614                 return map;
2615         res = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[0]));
2616         for (i = 1 ; i < map->n; ++i) {
2617                 struct isl_map *r2;
2618                 r2 = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[i]));
2619                 if (F_ISSET(map, ISL_MAP_DISJOINT))
2620                         res = isl_map_union_disjoint(res, r2);
2621                 else
2622                         res = isl_map_union(res, r2);
2623         }
2624         isl_map_free(map);
2625
2626         return res;
2627 }
2628
2629 struct isl_set *isl_basic_set_compute_divs(struct isl_basic_set *bset)
2630 {
2631         return (struct isl_set *)
2632                 isl_basic_map_compute_divs((struct isl_basic_map *)bset);
2633 }
2634
2635 struct isl_set *isl_set_compute_divs(struct isl_set *set)
2636 {
2637         return (struct isl_set *)
2638                 isl_map_compute_divs((struct isl_map *)set);
2639 }
2640
2641 struct isl_set *isl_map_domain(struct isl_map *map)
2642 {
2643         int i;
2644         struct isl_set *set;
2645
2646         if (!map)
2647                 goto error;
2648
2649         map = isl_map_cow(map);
2650         if (!map)
2651                 return NULL;
2652
2653         set = (struct isl_set *)map;
2654         set->dim = map->n_in;
2655         set->zero = 0;
2656         for (i = 0; i < map->n; ++i) {
2657                 set->p[i] = isl_basic_map_domain(map->p[i]);
2658                 if (!set->p[i])
2659                         goto error;
2660         }
2661         F_CLR(set, ISL_MAP_DISJOINT);
2662         return set;
2663 error:
2664         isl_map_free(map);
2665         return NULL;
2666 }
2667
2668 struct isl_map *isl_map_union_disjoint(
2669                         struct isl_map *map1, struct isl_map *map2)
2670 {
2671         int i;
2672         unsigned flags = 0;
2673         struct isl_map *map = NULL;
2674
2675         if (!map1 || !map2)
2676                 goto error;
2677
2678         if (map1->n == 0) {
2679                 isl_map_free(map1);
2680                 return map2;
2681         }
2682         if (map2->n == 0) {
2683                 isl_map_free(map2);
2684                 return map1;
2685         }
2686
2687         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
2688         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
2689         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
2690
2691         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
2692             F_ISSET(map2, ISL_MAP_DISJOINT))
2693                 FL_SET(flags, ISL_MAP_DISJOINT);
2694
2695         map = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
2696                                 map1->n + map2->n, flags);
2697         if (!map)
2698                 goto error;
2699         for (i = 0; i < map1->n; ++i) {
2700                 map = isl_map_add(map,
2701                                   isl_basic_map_copy(map1->p[i]));
2702                 if (!map)
2703                         goto error;
2704         }
2705         for (i = 0; i < map2->n; ++i) {
2706                 map = isl_map_add(map,
2707                                   isl_basic_map_copy(map2->p[i]));
2708                 if (!map)
2709                         goto error;
2710         }
2711         isl_map_free(map1);
2712         isl_map_free(map2);
2713         return map;
2714 error:
2715         isl_map_free(map);
2716         isl_map_free(map1);
2717         isl_map_free(map2);
2718         return NULL;
2719 }
2720
2721 struct isl_map *isl_map_union(struct isl_map *map1, struct isl_map *map2)
2722 {
2723         map1 = isl_map_union_disjoint(map1, map2);
2724         if (!map1)
2725                 return NULL;
2726         if (map1->n > 1)
2727                 F_CLR(map1, ISL_MAP_DISJOINT);
2728         return map1;
2729 }
2730
2731 struct isl_set *isl_set_union_disjoint(
2732                         struct isl_set *set1, struct isl_set *set2)
2733 {
2734         return (struct isl_set *)
2735                 isl_map_union_disjoint(
2736                         (struct isl_map *)set1, (struct isl_map *)set2);
2737 }
2738
2739 struct isl_set *isl_set_union(struct isl_set *set1, struct isl_set *set2)
2740 {
2741         return (struct isl_set *)
2742                 isl_map_union((struct isl_map *)set1, (struct isl_map *)set2);
2743 }
2744
2745 struct isl_map *isl_map_intersect_range(
2746                 struct isl_map *map, struct isl_set *set)
2747 {
2748         unsigned flags = 0;
2749         struct isl_map *result;
2750         int i, j;
2751
2752         if (!map || !set)
2753                 goto error;
2754
2755         if (F_ISSET(map, ISL_MAP_DISJOINT) &&
2756             F_ISSET(set, ISL_MAP_DISJOINT))
2757                 FL_SET(flags, ISL_MAP_DISJOINT);
2758
2759         result = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
2760                                 map->n * set->n, flags);
2761         if (!result)
2762                 goto error;
2763         for (i = 0; i < map->n; ++i)
2764                 for (j = 0; j < set->n; ++j) {
2765                         result = isl_map_add(result,
2766                             isl_basic_map_intersect_range(
2767                                 isl_basic_map_copy(map->p[i]),
2768                                 isl_basic_set_copy(set->p[j])));
2769                         if (!result)
2770                                 goto error;
2771                 }
2772         isl_map_free(map);
2773         isl_set_free(set);
2774         return result;
2775 error:
2776         isl_map_free(map);
2777         isl_set_free(set);
2778         return NULL;
2779 }
2780
2781 struct isl_map *isl_map_intersect_domain(
2782                 struct isl_map *map, struct isl_set *set)
2783 {
2784         return isl_map_reverse(
2785                 isl_map_intersect_range(isl_map_reverse(map), set));
2786 }
2787
2788 struct isl_map *isl_map_apply_domain(
2789                 struct isl_map *map1, struct isl_map *map2)
2790 {
2791         if (!map1 || !map2)
2792                 goto error;
2793         map1 = isl_map_reverse(map1);
2794         map1 = isl_map_apply_range(map1, map2);
2795         return isl_map_reverse(map1);
2796 error:
2797         isl_map_free(map1);
2798         isl_map_free(map2);
2799         return NULL;
2800 }
2801
2802 struct isl_map *isl_map_apply_range(
2803                 struct isl_map *map1, struct isl_map *map2)
2804 {
2805         struct isl_map *result;
2806         int i, j;
2807
2808         if (!map1 || !map2)
2809                 goto error;
2810
2811         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
2812         isl_assert(map1->ctx, map1->n_out == map2->n_in, goto error);
2813
2814         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map2->n_out,
2815                                 map1->n * map2->n, 0);
2816         if (!result)
2817                 goto error;
2818         for (i = 0; i < map1->n; ++i)
2819                 for (j = 0; j < map2->n; ++j) {
2820                         result = isl_map_add(result,
2821                             isl_basic_map_apply_range(
2822                                 isl_basic_map_copy(map1->p[i]),
2823                                 isl_basic_map_copy(map2->p[j])));
2824                         if (!result)
2825                                 goto error;
2826                 }
2827         isl_map_free(map1);
2828         isl_map_free(map2);
2829         if (result && result->n <= 1)
2830                 F_SET(result, ISL_MAP_DISJOINT);
2831         return result;
2832 error:
2833         isl_map_free(map1);
2834         isl_map_free(map2);
2835         return NULL;
2836 }
2837
2838 /*
2839  * returns range - domain
2840  */
2841 struct isl_basic_set *isl_basic_map_deltas(struct isl_basic_map *bmap)
2842 {
2843         struct isl_basic_set *bset;
2844         unsigned dim;
2845         int i;
2846
2847         if (!bmap)
2848                 goto error;
2849         isl_assert(bmap->ctx, bmap->n_in == bmap->n_out, goto error);
2850         dim = bmap->n_in;
2851         bset = isl_basic_set_from_basic_map(bmap);
2852         bset = isl_basic_set_extend(bset, bset->nparam, 3*dim, 0,
2853                                         dim, 0);
2854         bset = isl_basic_set_swap_vars(bset, 2*dim);
2855         for (i = 0; i < dim; ++i) {
2856                 int j = isl_basic_map_alloc_equality(
2857                                             (struct isl_basic_map *)bset);
2858                 if (j < 0)
2859                         goto error;
2860                 isl_seq_clr(bset->eq[j],
2861                             1 + bset->nparam + bset->dim + bset->extra);
2862                 isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
2863                 isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
2864                 isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
2865         }
2866         return isl_basic_set_project_out(bset, 2*dim, 0);
2867 error:
2868         isl_basic_map_free(bmap);
2869         return NULL;
2870 }
2871
2872 /* If the only constraints a div d=floor(f/m)
2873  * appears in are its two defining constraints
2874  *
2875  *      f - m d >=0
2876  *      -(f - (m - 1)) + m d >= 0
2877  *
2878  * then it can safely be removed.
2879  */
2880 static int div_is_redundant(struct isl_basic_map *bmap, int div)
2881 {
2882         int i;
2883         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2884
2885         for (i = 0; i < bmap->n_eq; ++i)
2886                 if (!isl_int_is_zero(bmap->eq[i][pos]))
2887                         return 0;
2888
2889         for (i = 0; i < bmap->n_ineq; ++i) {
2890                 if (isl_int_is_zero(bmap->ineq[i][pos]))
2891                         continue;
2892                 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2893                         int neg;
2894                         isl_int_sub(bmap->div[div][1],
2895                                         bmap->div[div][1], bmap->div[div][0]);
2896                         isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
2897                         neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
2898                         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
2899                         isl_int_add(bmap->div[div][1],
2900                                         bmap->div[div][1], bmap->div[div][0]);
2901                         if (!neg)
2902                                 return 0;
2903                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2904                                                     bmap->n_div-div-1) != -1)
2905                                 return 0;
2906                 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2907                         if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
2908                                 return 0;
2909                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2910                                                     bmap->n_div-div-1) != -1)
2911                                 return 0;
2912                 } else
2913                         return 0;
2914         }
2915
2916         for (i = 0; i < bmap->n_div; ++i)
2917                 if (!isl_int_is_zero(bmap->div[i][1+pos]))
2918                         return 0;
2919
2920         return 1;
2921 }
2922
2923 /*
2924  * Remove divs that don't occur in any of the constraints or other divs.
2925  * These can arise when dropping some of the variables in a quast
2926  * returned by piplib.
2927  */
2928 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
2929 {
2930         int i;
2931
2932         if (!bmap)
2933                 return NULL;
2934
2935         for (i = bmap->n_div-1; i >= 0; --i) {
2936                 if (!div_is_redundant(bmap, i))
2937                         continue;
2938                 bmap = isl_basic_map_drop_div(bmap, i);
2939         }
2940         return bmap;
2941 }
2942
2943 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
2944 {
2945         bmap = remove_redundant_divs(bmap);
2946         if (!bmap)
2947                 return NULL;
2948         F_SET(bmap, ISL_BASIC_SET_FINAL);
2949         return bmap;
2950 }
2951
2952 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
2953 {
2954         return (struct isl_basic_set *)
2955                 isl_basic_map_finalize((struct isl_basic_map *)bset);
2956 }
2957
2958 struct isl_set *isl_set_finalize(struct isl_set *set)
2959 {
2960         int i;
2961
2962         if (!set)
2963                 return NULL;
2964         for (i = 0; i < set->n; ++i) {
2965                 set->p[i] = isl_basic_set_finalize(set->p[i]);
2966                 if (!set->p[i])
2967                         goto error;
2968         }
2969         return set;
2970 error:
2971         isl_set_free(set);
2972         return NULL;
2973 }
2974
2975 struct isl_map *isl_map_finalize(struct isl_map *map)
2976 {
2977         int i;
2978
2979         if (!map)
2980                 return NULL;
2981         for (i = 0; i < map->n; ++i) {
2982                 map->p[i] = isl_basic_map_finalize(map->p[i]);
2983                 if (!map->p[i])
2984                         goto error;
2985         }
2986         return map;
2987 error:
2988         isl_map_free(map);
2989         return NULL;
2990 }
2991
2992 struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
2993                 unsigned nparam, unsigned dim)
2994 {
2995         struct isl_basic_map *bmap;
2996         int i;
2997
2998         bmap = isl_basic_map_alloc(ctx, nparam, dim, dim, 0, dim, 0);
2999         if (!bmap)
3000                 goto error;
3001
3002         for (i = 0; i < dim; ++i) {
3003                 int j = isl_basic_map_alloc_equality(bmap);
3004                 if (j < 0)
3005                         goto error;
3006                 isl_seq_clr(bmap->eq[j],
3007                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
3008                 isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
3009                 isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
3010         }
3011         return isl_basic_map_finalize(bmap);
3012 error:
3013         isl_basic_map_free(bmap);
3014         return NULL;
3015 }
3016
3017 struct isl_map *isl_map_identity(struct isl_ctx *ctx,
3018                 unsigned nparam, unsigned dim)
3019 {
3020         struct isl_map *map = isl_map_alloc(ctx, nparam, dim, dim, 1,
3021                                                 ISL_MAP_DISJOINT);
3022         if (!map)
3023                 goto error;
3024         map = isl_map_add(map,
3025                 isl_basic_map_identity(ctx, nparam, dim));
3026         return map;
3027 error:
3028         isl_map_free(map);
3029         return NULL;
3030 }
3031
3032 int isl_set_is_equal(struct isl_set *set1, struct isl_set *set2)
3033 {
3034         return isl_map_is_equal((struct isl_map *)set1, (struct isl_map *)set2);
3035 }
3036
3037 int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
3038 {
3039         return isl_map_is_subset(
3040                         (struct isl_map *)set1, (struct isl_map *)set2);
3041 }
3042
3043 int isl_basic_map_is_subset(
3044                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3045 {
3046         int is_subset;
3047         struct isl_map *map1;
3048         struct isl_map *map2;
3049
3050         if (!bmap1 || !bmap2)
3051                 return -1;
3052
3053         map1 = isl_map_from_basic_map(isl_basic_map_copy(bmap1));
3054         map2 = isl_map_from_basic_map(isl_basic_map_copy(bmap2));
3055
3056         is_subset = isl_map_is_subset(map1, map2);
3057
3058         isl_map_free(map1);
3059         isl_map_free(map2);
3060
3061         return is_subset;
3062 }
3063
3064 int isl_basic_map_is_equal(
3065                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3066 {
3067         int is_subset;
3068
3069         if (!bmap1 || !bmap2)
3070                 return -1;
3071         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3072         if (is_subset != 1)
3073                 return is_subset;
3074         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3075         return is_subset;
3076 }
3077
3078 int isl_basic_set_is_equal(
3079                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3080 {
3081         return isl_basic_map_is_equal(
3082                 (struct isl_basic_map *)bset1, (struct isl_basic_map *)bset2);
3083 }
3084
3085 int isl_map_is_empty(struct isl_map *map)
3086 {
3087         int i;
3088         int is_empty;
3089
3090         if (!map)
3091                 return -1;
3092         for (i = 0; i < map->n; ++i) {
3093                 is_empty = isl_basic_map_is_empty(map->p[i]);
3094                 if (is_empty < 0)
3095                         return -1;
3096                 if (!is_empty)
3097                         return 0;
3098         }
3099         return 1;
3100 }
3101
3102 int isl_set_is_empty(struct isl_set *set)
3103 {
3104         return isl_map_is_empty((struct isl_map *)set);
3105 }
3106
3107 int isl_map_is_subset(struct isl_map *map1, struct isl_map *map2)
3108 {
3109         int i;
3110         int is_subset = 0;
3111         struct isl_map *diff;
3112
3113         if (!map1 || !map2)
3114                 return -1;
3115
3116         if (isl_map_is_empty(map1))
3117                 return 1;
3118
3119         if (isl_map_is_empty(map2))
3120                 return 0;
3121
3122         diff = isl_map_subtract(isl_map_copy(map1), isl_map_copy(map2));
3123         if (!diff)
3124                 return -1;
3125
3126         is_subset = isl_map_is_empty(diff);
3127         isl_map_free(diff);
3128
3129         return is_subset;
3130 }
3131
3132 int isl_map_is_equal(struct isl_map *map1, struct isl_map *map2)
3133 {
3134         int is_subset;
3135
3136         if (!map1 || !map2)
3137                 return -1;
3138         is_subset = isl_map_is_subset(map1, map2);
3139         if (is_subset != 1)
3140                 return is_subset;
3141         is_subset = isl_map_is_subset(map2, map1);
3142         return is_subset;
3143 }
3144
3145 int isl_basic_map_is_strict_subset(
3146                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3147 {
3148         int is_subset;
3149
3150         if (!bmap1 || !bmap2)
3151                 return -1;
3152         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3153         if (is_subset != 1)
3154                 return is_subset;
3155         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3156         if (is_subset == -1)
3157                 return is_subset;
3158         return !is_subset;
3159 }
3160
3161 static int basic_map_contains(struct isl_basic_map *bmap, struct isl_vec *vec)
3162 {
3163         int i;
3164         unsigned total;
3165         isl_int s;
3166
3167         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3168         if (total != vec->size)
3169                 return -1;
3170
3171         isl_int_init(s);
3172
3173         for (i = 0; i < bmap->n_eq; ++i) {
3174                 isl_seq_inner_product(vec->block.data, bmap->eq[i], total, &s);
3175                 if (!isl_int_is_zero(s)) {
3176                         isl_int_clear(s);
3177                         return 0;
3178                 }
3179         }
3180
3181         for (i = 0; i < bmap->n_ineq; ++i) {
3182                 isl_seq_inner_product(vec->block.data, bmap->ineq[i], total, &s);
3183                 if (isl_int_is_neg(s)) {
3184                         isl_int_clear(s);
3185                         return 0;
3186                 }
3187         }
3188
3189         isl_int_clear(s);
3190
3191         return 1;
3192 }
3193
3194 int isl_basic_map_is_universe(struct isl_basic_map *bmap)
3195 {
3196         if (!bmap)
3197                 return -1;
3198         return bmap->n_eq == 0 && bmap->n_ineq == 0;
3199 }
3200
3201 int isl_basic_map_is_empty(struct isl_basic_map *bmap)
3202 {
3203         struct isl_basic_set *bset = NULL;
3204         struct isl_vec *sample = NULL;
3205         int empty;
3206         unsigned total;
3207
3208         if (!bmap)
3209                 return -1;
3210
3211         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
3212                 return 1;
3213
3214         if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
3215                 struct isl_basic_map *copy = isl_basic_map_copy(bmap);
3216                 copy = isl_basic_map_convex_hull(copy);
3217                 empty = F_ISSET(copy, ISL_BASIC_MAP_EMPTY);
3218                 isl_basic_map_free(copy);
3219                 return empty;
3220         }
3221
3222         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3223         if (bmap->sample && bmap->sample->size == total) {
3224                 int contains = basic_map_contains(bmap, bmap->sample);
3225                 if (contains < 0)
3226                         return -1;
3227                 if (contains)
3228                         return 0;
3229         }
3230         bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
3231         if (!bset)
3232                 return -1;
3233         sample = isl_basic_set_sample(bset);
3234         if (!sample)
3235                 return -1;
3236         empty = sample->size == 0;
3237         if (bmap->sample)
3238                 isl_vec_free(bmap->ctx, bmap->sample);
3239         bmap->sample = sample;
3240
3241         return empty;
3242 }
3243
3244 struct isl_map *isl_basic_map_union(
3245         struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3246 {
3247         struct isl_map *map;
3248         if (!bmap1 || !bmap2)
3249                 return NULL;
3250
3251         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
3252         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
3253         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
3254
3255         map = isl_map_alloc(bmap1->ctx, bmap1->nparam,
3256                                 bmap1->n_in, bmap1->n_out, 2, 0);
3257         if (!map)
3258                 goto error;
3259         map = isl_map_add(map, bmap1);
3260         map = isl_map_add(map, bmap2);
3261         return map;
3262 error:
3263         isl_basic_map_free(bmap1);
3264         isl_basic_map_free(bmap2);
3265         return NULL;
3266 }
3267
3268 struct isl_set *isl_basic_set_union(
3269                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3270 {
3271         return (struct isl_set *)isl_basic_map_union(
3272                                             (struct isl_basic_map *)bset1,
3273                                             (struct isl_basic_map *)bset2);
3274 }
3275
3276 /* Order divs such that any div only depends on previous divs */
3277 static struct isl_basic_map *order_divs(struct isl_basic_map *bmap)
3278 {
3279         int i;
3280         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
3281
3282         for (i = 0; i < bmap->n_div; ++i) {
3283                 int pos;
3284                 pos = isl_seq_first_non_zero(bmap->div[i]+1+1+off+i,
3285                                                             bmap->n_div-i);
3286                 if (pos == -1)
3287                         continue;
3288                 swap_div(bmap, i, pos);
3289                 --i;
3290         }
3291         return bmap;
3292 }
3293
3294 static int find_div(struct isl_basic_map *dst,
3295                         struct isl_basic_map *src, unsigned div)
3296 {
3297         int i;
3298
3299         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3300
3301         for (i = 0; i < dst->n_div; ++i)
3302                 if (isl_seq_eq(dst->div[i], src->div[div], 1+1+total) &&
3303                     isl_seq_first_non_zero(dst->div[i]+1+1+total,
3304                                                 dst->n_div - src->n_div) == -1)
3305                         return i;
3306         return -1;
3307 }
3308
3309 struct isl_basic_map *isl_basic_map_align_divs(
3310                 struct isl_basic_map *dst, struct isl_basic_map *src)
3311 {
3312         int i;
3313         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3314
3315         if (src->n_div == 0)
3316                 return dst;
3317
3318         src = order_divs(src);
3319         dst = isl_basic_map_extend(dst, dst->nparam, dst->n_in,
3320                         dst->n_out, src->n_div, 0, 2 * src->n_div);
3321         if (!dst)
3322                 return NULL;
3323         for (i = 0; i < src->n_div; ++i) {
3324                 int j = find_div(dst, src, i);
3325                 if (j < 0) {
3326                         j = isl_basic_map_alloc_div(dst);
3327                         if (j < 0)
3328                                 goto error;
3329                         isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
3330                         isl_seq_clr(dst->div[j]+1+1+total,
3331                                                     dst->extra - src->n_div);
3332                         if (add_div_constraints(dst, j) < 0)
3333                                 goto error;
3334                 }
3335                 if (j != i)
3336                         swap_div(dst, i, j);
3337         }
3338         return dst;
3339 error:
3340         isl_basic_map_free(dst);
3341         return NULL;
3342 }
3343
3344 struct isl_map *isl_map_align_divs(struct isl_map *map)
3345 {
3346         int i;
3347
3348         map = isl_map_compute_divs(map);
3349         map = isl_map_cow(map);
3350         if (!map)
3351                 return NULL;
3352
3353         for (i = 1; i < map->n; ++i)
3354                 map->p[0] = isl_basic_map_align_divs(map->p[0], map->p[i]);
3355         for (i = 1; i < map->n; ++i)
3356                 map->p[i] = isl_basic_map_align_divs(map->p[i], map->p[0]);
3357
3358         return map;
3359 }
3360
3361 static struct isl_map *add_cut_constraint(struct isl_map *dst,
3362                 struct isl_basic_map *src, isl_int *c,
3363                 unsigned len, int oppose)
3364 {
3365         struct isl_basic_map *copy = NULL;
3366         int is_empty;
3367         int k;
3368         unsigned total;
3369
3370         copy = isl_basic_map_copy(src);
3371         copy = isl_basic_map_cow(copy);
3372         if (!copy)
3373                 goto error;
3374         copy = isl_basic_map_extend(copy,
3375                 copy->nparam, copy->n_in, copy->n_out, 0, 0, 1);
3376         k = isl_basic_map_alloc_inequality(copy);
3377         if (k < 0)
3378                 goto error;
3379         if (oppose)
3380                 isl_seq_neg(copy->ineq[k], c, len);
3381         else
3382                 isl_seq_cpy(copy->ineq[k], c, len);
3383         total = 1 + copy->nparam + copy->n_in + copy->n_out + copy->extra;
3384         isl_seq_clr(copy->ineq[k]+len, total - len);
3385         isl_inequality_negate(copy, k);
3386         copy = isl_basic_map_simplify(copy);
3387         copy = isl_basic_map_finalize(copy);
3388         is_empty = isl_basic_map_is_empty(copy);
3389         if (is_empty < 0)
3390                 goto error;
3391         if (!is_empty)
3392                 dst = isl_map_add(dst, copy);
3393         else
3394                 isl_basic_map_free(copy);
3395         return dst;
3396 error:
3397         isl_basic_map_free(copy);
3398         isl_map_free(dst);
3399         return NULL;
3400 }
3401
3402 static struct isl_map *subtract(struct isl_map *map, struct isl_basic_map *bmap)
3403 {
3404         int i, j, k;
3405         unsigned flags = 0;
3406         struct isl_map *rest = NULL;
3407         unsigned max;
3408         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3409
3410         assert(bmap);
3411
3412         if (!map)
3413                 goto error;
3414
3415         if (F_ISSET(map, ISL_MAP_DISJOINT))
3416                 FL_SET(flags, ISL_MAP_DISJOINT);
3417
3418         max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
3419         rest = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
3420                                 max, flags);
3421         if (!rest)
3422                 goto error;
3423
3424         for (i = 0; i < map->n; ++i) {
3425                 map->p[i] = isl_basic_map_align_divs(map->p[i], bmap);
3426                 if (!map->p[i])
3427                         goto error;
3428         }
3429
3430         for (j = 0; j < map->n; ++j)
3431                 map->p[j] = isl_basic_map_cow(map->p[j]);
3432
3433         for (i = 0; i < bmap->n_eq; ++i) {
3434                 for (j = 0; j < map->n; ++j) {
3435                         rest = add_cut_constraint(rest,
3436                                 map->p[j], bmap->eq[i], 1+total, 0);
3437                         if (!rest)
3438                                 goto error;
3439
3440                         rest = add_cut_constraint(rest,
3441                                 map->p[j], bmap->eq[i], 1+total, 1);
3442                         if (!rest)
3443                                 goto error;
3444
3445                         map->p[j] = isl_basic_map_extend(map->p[j],
3446                                 map->p[j]->nparam, map->p[j]->n_in,
3447                                 map->p[j]->n_out, 0, 1, 0);
3448                         if (!map->p[j])
3449                                 goto error;
3450                         k = isl_basic_map_alloc_equality(map->p[j]);
3451                         if (k < 0)
3452                                 goto error;
3453                         isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
3454                         isl_seq_clr(map->p[j]->eq[k]+1+total,
3455                                         map->p[j]->extra - bmap->n_div);
3456                 }
3457         }
3458
3459         for (i = 0; i < bmap->n_ineq; ++i) {
3460                 for (j = 0; j < map->n; ++j) {
3461                         rest = add_cut_constraint(rest,
3462                                 map->p[j], bmap->ineq[i], 1+total, 0);
3463                         if (!rest)
3464                                 goto error;
3465
3466                         map->p[j] = isl_basic_map_extend(map->p[j],
3467                                 map->p[j]->nparam, map->p[j]->n_in,
3468                                 map->p[j]->n_out, 0, 0, 1);
3469                         if (!map->p[j])
3470                                 goto error;
3471                         k = isl_basic_map_alloc_inequality(map->p[j]);
3472                         if (k < 0)
3473                                 goto error;
3474                         isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
3475                         isl_seq_clr(map->p[j]->ineq[k]+1+total,
3476                                         map->p[j]->extra - bmap->n_div);
3477                 }
3478         }
3479
3480         isl_map_free(map);
3481         return rest;
3482 error:
3483         isl_map_free(map);
3484         isl_map_free(rest);
3485         return NULL;
3486 }
3487
3488 struct isl_map *isl_map_subtract(struct isl_map *map1, struct isl_map *map2)
3489 {
3490         int i;
3491         if (!map1 || !map2)
3492                 goto error;
3493
3494         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
3495         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
3496         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
3497
3498         if (isl_map_is_empty(map2)) {
3499                 isl_map_free(map2);
3500                 return map1;
3501         }
3502
3503         map1 = isl_map_compute_divs(map1);
3504         map2 = isl_map_compute_divs(map2);
3505         if (!map1 || !map2)
3506                 goto error;
3507
3508         for (i = 0; map1 && i < map2->n; ++i)
3509                 map1 = subtract(map1, map2->p[i]);
3510
3511         isl_map_free(map2);
3512         return map1;
3513 error:
3514         isl_map_free(map1);
3515         isl_map_free(map2);
3516         return NULL;
3517 }
3518
3519 struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
3520 {
3521         return (struct isl_set *)
3522                 isl_map_subtract(
3523                         (struct isl_map *)set1, (struct isl_map *)set2);
3524 }
3525
3526 struct isl_set *isl_set_apply(struct isl_set *set, struct isl_map *map)
3527 {
3528         if (!set || !map)
3529                 goto error;
3530         isl_assert(set->ctx, set->dim == map->n_in, goto error);
3531         map = isl_map_intersect_domain(map, set);
3532         set = isl_map_range(map);
3533         return set;
3534 error:
3535         isl_set_free(set);
3536         isl_map_free(map);
3537         return NULL;
3538 }
3539
3540 /* There is no need to cow as removing empty parts doesn't change
3541  * the meaning of the set.
3542  */
3543 struct isl_map *isl_map_remove_empty_parts(struct isl_map *map)
3544 {
3545         int i;
3546
3547         if (!map)
3548                 return NULL;
3549
3550         for (i = map->n-1; i >= 0; --i) {
3551                 if (!F_ISSET(map->p[i], ISL_BASIC_MAP_EMPTY))
3552                         continue;
3553                 isl_basic_map_free(map->p[i]);
3554                 if (i != map->n-1)
3555                         map->p[i] = map->p[map->n-1];
3556                 map->n--;
3557         }
3558
3559         return map;
3560 }
3561
3562 struct isl_set *isl_set_remove_empty_parts(struct isl_set *set)
3563 {
3564         return (struct isl_set *)
3565                 isl_map_remove_empty_parts((struct isl_map *)set);
3566 }
3567
3568 struct isl_basic_set *isl_set_copy_basic_set(struct isl_set *set)
3569 {
3570         struct isl_basic_set *bset;
3571         if (!set || set->n == 0)
3572                 return NULL;
3573         bset = set->p[set->n-1];
3574         isl_assert(set->ctx, F_ISSET(bset, ISL_BASIC_SET_FINAL), return NULL);
3575         return isl_basic_set_copy(bset);
3576 }
3577
3578 struct isl_set *isl_set_drop_basic_set(struct isl_set *set,
3579                                                 struct isl_basic_set *bset)
3580 {
3581         int i;
3582
3583         if (!set || !bset)
3584                 goto error;
3585         for (i = set->n-1; i >= 0; --i) {
3586                 if (set->p[i] != bset)
3587                         continue;
3588                 set = isl_set_cow(set);
3589                 if (!set)
3590                         goto error;
3591                 isl_basic_set_free(set->p[i]);
3592                 if (i != set->n-1)
3593                         set->p[i] = set->p[set->n-1];
3594                 set->n--;
3595                 return set;
3596         }
3597         isl_basic_set_free(bset);
3598         return set;
3599 error:
3600         isl_set_free(set);
3601         isl_basic_set_free(bset);
3602         return NULL;
3603 }
3604
3605 /* Given two _disjoint_ basic sets bset1 and bset2, check whether
3606  * for any common value of the parameters and dimensions preceding dim
3607  * in both basic sets, the values of dimension pos in bset1 are
3608  * smaller or larger then those in bset2.
3609  *
3610  * Returns
3611  *       1 if bset1 follows bset2
3612  *      -1 if bset1 precedes bset2
3613  *       0 if bset1 and bset2 are incomparable
3614  *      -2 if some error occurred.
3615  */
3616 int isl_basic_set_compare_at(struct isl_basic_set *bset1,
3617         struct isl_basic_set *bset2, int pos)
3618 {
3619         struct isl_basic_map *bmap1 = NULL;
3620         struct isl_basic_map *bmap2 = NULL;
3621         struct isl_ctx *ctx;
3622         struct isl_vec *obj;
3623         unsigned total;
3624         isl_int num, den;
3625         enum isl_lp_result res;
3626         int cmp;
3627
3628         if (!bset1 || !bset2)
3629                 return -2;
3630
3631         bmap1 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset1),
3632                                                 pos, bset1->dim - pos);
3633         bmap2 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset2),
3634                                                 pos, bset2->dim - pos);
3635         if (!bmap1 || !bmap2)
3636                 goto error;
3637         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
3638                         bmap1->n_in, bmap1->n_out + bmap2->n_out,
3639                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
3640         if (!bmap1)
3641                 goto error;
3642         total = bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
3643         ctx = bmap1->ctx;
3644         obj = isl_vec_alloc(ctx, total);
3645         isl_seq_clr(obj->block.data, total);
3646         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in], 1);
3647         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in+
3648                                         bmap1->n_out-bmap2->n_out], -1);
3649         if (!obj)
3650                 goto error;
3651         bmap1 = add_constraints(bmap1, bmap2, 0, bmap1->n_out - bmap2->n_out);
3652         isl_int_init(num);
3653         isl_int_init(den);
3654         res = isl_solve_lp(bmap1, 0, obj->block.data, ctx->one, &num, &den);
3655         if (res == isl_lp_empty)
3656                 cmp = 0;
3657         else if (res == isl_lp_ok && isl_int_is_pos(num))
3658                 cmp = 1;
3659         else if ((res == isl_lp_ok && isl_int_is_neg(num)) ||
3660                   res == isl_lp_unbounded)
3661                 cmp = -1;
3662         else
3663                 cmp = -2;
3664         isl_int_clear(num);
3665         isl_int_clear(den);
3666         isl_basic_map_free(bmap1);
3667         isl_vec_free(ctx, obj);
3668         return cmp;
3669 error:
3670         isl_basic_map_free(bmap1);
3671         isl_basic_map_free(bmap2);
3672         return -2;
3673 }