isl_map.c

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