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 struct isl_basic_map *remove_duplicate_divs(
1111         struct isl_basic_map *bmap, int *progress)
1112 {
1113         unsigned int size;
1114         int *index;
1115         int k, l, h;
1116         int bits;
1117         struct isl_blk eq;
1118         unsigned total_var = bmap->nparam + bmap->n_in + bmap->n_out;
1119         unsigned total = total_var + bmap->n_div;
1120         struct isl_ctx *ctx;
1121
1122         if (bmap->n_div <= 1)
1123                 return bmap;
1124
1125         ctx = bmap->ctx;
1126         for (k = bmap->n_div - 1; k >= 0; --k)
1127                 if (!isl_int_is_zero(bmap->div[k][0]))
1128                         break;
1129         if (k <= 0)
1130                 return bmap;
1131
1132         size = round_up(4 * bmap->n_div / 3 - 1);
1133         bits = ffs(size) - 1;
1134         index = isl_alloc_array(ctx, int, size);
1135         memset(index, 0, size * sizeof(int));
1136         if (!index)
1137                 return bmap;
1138         eq = isl_blk_alloc(ctx, 1+total);
1139         if (isl_blk_is_error(eq))
1140                 goto out;
1141
1142         isl_seq_clr(eq.data, 1+total);
1143         index[isl_seq_hash(bmap->div[k], 2+total, bits)] = k + 1;
1144         for (--k; k >= 0; --k) {
1145                 u_int32_t hash;
1146
1147                 if (isl_int_is_zero(bmap->div[k][0]))
1148                         continue;
1149
1150                 hash = isl_seq_hash(bmap->div[k], 2+total, bits);
1151                 for (h = hash; index[h]; h = (h+1) % size)
1152                         if (isl_seq_eq(bmap->div[k],
1153                                        bmap->div[index[h]-1], 2+total))
1154                                 break;
1155                 if (index[h]) {
1156                         *progress = 1;
1157                         l = index[h] - 1;
1158                         isl_int_set_si(eq.data[1+total_var+k], -1);
1159                         isl_int_set_si(eq.data[1+total_var+l], 1);
1160                         eliminate_div(bmap, eq.data, l);
1161                         isl_int_set_si(eq.data[1+total_var+k], 0);
1162                         isl_int_set_si(eq.data[1+total_var+l], 0);
1163                 }
1164                 index[h] = k+1;
1165         }
1166
1167         isl_blk_free(ctx, eq);
1168 out:
1169         free(index);
1170         return bmap;
1171 }
1172
1173 /* Elimininate divs based on equalities
1174  */
1175 static struct isl_basic_map *eliminate_divs_eq(
1176                 struct isl_basic_map *bmap, int *progress)
1177 {
1178         int d;
1179         int i;
1180         int modified = 0;
1181         unsigned off;
1182
1183         if (!bmap)
1184                 return NULL;
1185
1186         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1187
1188         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1189                 for (i = 0; i < bmap->n_eq; ++i) {
1190                         if (!isl_int_is_one(bmap->eq[i][off + d]) &&
1191                             !isl_int_is_negone(bmap->eq[i][off + d]))
1192                                 continue;
1193                         modified = 1;
1194                         *progress = 1;
1195                         eliminate_div(bmap, bmap->eq[i], d);
1196                         isl_basic_map_drop_equality(bmap, i);
1197                         break;
1198                 }
1199         }
1200         if (modified)
1201                 return eliminate_divs_eq(bmap, progress);
1202         return bmap;
1203 }
1204
1205 static struct isl_basic_map *normalize_constraints(struct isl_basic_map *bmap)
1206 {
1207         int i;
1208         isl_int gcd;
1209         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1210
1211         isl_int_init(gcd);
1212         for (i = bmap->n_eq - 1; i >= 0; --i) {
1213                 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
1214                 if (isl_int_is_zero(gcd)) {
1215                         if (!isl_int_is_zero(bmap->eq[i][0])) {
1216                                 bmap = isl_basic_map_set_to_empty(bmap);
1217                                 break;
1218                         }
1219                         isl_basic_map_drop_equality(bmap, i);
1220                         continue;
1221                 }
1222                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1223                         isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
1224                 if (isl_int_is_one(gcd))
1225                         continue;
1226                 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
1227                         bmap = isl_basic_map_set_to_empty(bmap);
1228                         break;
1229                 }
1230                 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
1231         }
1232
1233         for (i = bmap->n_ineq - 1; i >= 0; --i) {
1234                 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
1235                 if (isl_int_is_zero(gcd)) {
1236                         if (isl_int_is_neg(bmap->ineq[i][0])) {
1237                                 bmap = isl_basic_map_set_to_empty(bmap);
1238                                 break;
1239                         }
1240                         isl_basic_map_drop_inequality(bmap, i);
1241                         continue;
1242                 }
1243                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1244                         isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
1245                 if (isl_int_is_one(gcd))
1246                         continue;
1247                 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
1248                 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
1249         }
1250         isl_int_clear(gcd);
1251
1252         return bmap;
1253 }
1254
1255
1256 /* Eliminate the specified variables from the constraints using
1257  * Fourier-Motzkin.  The variables themselves are not removed.
1258  */
1259 struct isl_basic_map *isl_basic_map_eliminate_vars(
1260         struct isl_basic_map *bmap, int pos, unsigned n)
1261 {
1262         int d;
1263         int i, j, k;
1264         unsigned total;
1265         unsigned extra;
1266
1267         if (!bmap)
1268                 return NULL;
1269         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1270         extra = bmap->extra - bmap->n_div;
1271
1272         bmap = isl_basic_map_cow(bmap);
1273         for (d = pos + n - 1; d >= pos; --d) {
1274                 int n_lower, n_upper;
1275                 if (!bmap)
1276                         return NULL;
1277                 for (i = 0; i < bmap->n_eq; ++i) {
1278                         if (isl_int_is_zero(bmap->eq[i][1+d]))
1279                                 continue;
1280                         eliminate_var_using_equality(bmap, d, bmap->eq[i], NULL);
1281                         isl_basic_map_drop_equality(bmap, i);
1282                         break;
1283                 }
1284                 if (i < bmap->n_eq)
1285                         continue;
1286                 n_lower = 0;
1287                 n_upper = 0;
1288                 for (i = 0; i < bmap->n_ineq; ++i) {
1289                         if (isl_int_is_pos(bmap->ineq[i][1+d]))
1290                                 n_lower++;
1291                         else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1292                                 n_upper++;
1293                 }
1294                 bmap = isl_basic_map_extend(bmap,
1295                                 bmap->nparam, bmap->n_in, bmap->n_out, 0,
1296                                 0, n_lower * n_upper);
1297                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1298                         int last;
1299                         if (isl_int_is_zero(bmap->ineq[i][1+d]))
1300                                 continue;
1301                         last = -1;
1302                         for (j = 0; j < i; ++j) {
1303                                 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1304                                         continue;
1305                                 last = j;
1306                                 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1307                                     isl_int_sgn(bmap->ineq[j][1+d]))
1308                                         continue;
1309                                 k = isl_basic_map_alloc_inequality(bmap);
1310                                 if (k < 0)
1311                                         goto error;
1312                                 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1313                                                 1+total);
1314                                 isl_seq_clr(bmap->ineq[k]+1+total, extra);
1315                                 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1316                                                 1+d, 1+total, NULL);
1317                         }
1318                         isl_basic_map_drop_inequality(bmap, i);
1319                         i = last + 1;
1320                 }
1321                 if (n_lower > 0 && n_upper > 0) {
1322                         bmap = normalize_constraints(bmap);
1323                         bmap = remove_duplicate_constraints(bmap, NULL);
1324                         bmap = isl_basic_map_gauss(bmap, NULL);
1325                         bmap = isl_basic_map_convex_hull(bmap);
1326                         if (!bmap)
1327                                 goto error;
1328                         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1329                                 break;
1330                 }
1331         }
1332         F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1333         return bmap;
1334 error:
1335         isl_basic_map_free(bmap);
1336         return NULL;
1337 }
1338
1339 struct isl_basic_set *isl_basic_set_eliminate_vars(
1340         struct isl_basic_set *bset, unsigned pos, unsigned n)
1341 {
1342         return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1343                         (struct isl_basic_map *)bset, pos, n);
1344 }
1345
1346 /* Eliminate the specified n dimensions starting at first from the
1347  * constraints using Fourier-Motzkin, The dimensions themselves
1348  * are not removed.
1349  */
1350 struct isl_set *isl_set_eliminate_dims(struct isl_set *set,
1351         unsigned first, unsigned n)
1352 {
1353         int i;
1354
1355         if (n == 0)
1356                 return set;
1357
1358         set = isl_set_cow(set);
1359         if (!set)
1360                 return NULL;
1361         isl_assert(set->ctx, first+n <= set->dim, goto error);
1362
1363         for (i = 0; i < set->n; ++i) {
1364                 set->p[i] = isl_basic_set_eliminate_vars(set->p[i],
1365                                                         set->nparam + first, n);
1366                 if (!set->p[i])
1367                         goto error;
1368         }
1369         return set;
1370 error:
1371         isl_set_free(set);
1372         return NULL;
1373 }
1374
1375 /* Project out n dimensions starting at first using Fourier-Motzkin */
1376 struct isl_set *isl_set_remove_dims(struct isl_set *set,
1377         unsigned first, unsigned n)
1378 {
1379         set = isl_set_eliminate_dims(set, first, n);
1380         set = isl_set_drop_dims(set, first, n);
1381         return set;
1382 }
1383
1384 /* Project out n inputs starting at first using Fourier-Motzkin */
1385 struct isl_map *isl_map_remove_inputs(struct isl_map *map,
1386         unsigned first, unsigned n)
1387 {
1388         int i;
1389
1390         if (n == 0)
1391                 return map;
1392
1393         map = isl_map_cow(map);
1394         if (!map)
1395                 return NULL;
1396         isl_assert(map->ctx, first+n <= map->n_in, goto error);
1397
1398         for (i = 0; i < map->n; ++i) {
1399                 map->p[i] = isl_basic_map_eliminate_vars(map->p[i],
1400                                                         map->nparam + first, n);
1401                 if (!map->p[i])
1402                         goto error;
1403         }
1404         map = isl_map_drop_inputs(map, first, n);
1405         return map;
1406 error:
1407         isl_map_free(map);
1408         return NULL;
1409 }
1410
1411 /* Project out n dimensions starting at first using Fourier-Motzkin */
1412 struct isl_basic_set *isl_basic_set_remove_dims(struct isl_basic_set *bset,
1413         unsigned first, unsigned n)
1414 {
1415         bset = isl_basic_set_eliminate_vars(bset, bset->nparam + first, n);
1416         bset = isl_basic_set_drop_dims(bset, first, n);
1417         return bset;
1418 }
1419
1420 /* Elimininate divs based on inequalities
1421  */
1422 static struct isl_basic_map *eliminate_divs_ineq(
1423                 struct isl_basic_map *bmap, int *progress)
1424 {
1425         int d;
1426         int i;
1427         unsigned off;
1428         struct isl_ctx *ctx;
1429
1430         if (!bmap)
1431                 return NULL;
1432
1433         ctx = bmap->ctx;
1434         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1435
1436         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1437                 for (i = 0; i < bmap->n_eq; ++i)
1438                         if (!isl_int_is_zero(bmap->eq[i][off + d]))
1439                                 break;
1440                 if (i < bmap->n_eq)
1441                         continue;
1442                 for (i = 0; i < bmap->n_ineq; ++i)
1443                         if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
1444                                 break;
1445                 if (i < bmap->n_ineq)
1446                         continue;
1447                 *progress = 1;
1448                 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
1449                 if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1450                         break;
1451                 bmap = isl_basic_map_drop_div(bmap, d);
1452                 if (!bmap)
1453                         break;
1454         }
1455         return bmap;
1456 }
1457
1458 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1459 {
1460         int progress = 1;
1461         if (!bmap)
1462                 return NULL;
1463         while (progress) {
1464                 progress = 0;
1465                 bmap = normalize_constraints(bmap);
1466                 bmap = eliminate_divs_eq(bmap, &progress);
1467                 bmap = eliminate_divs_ineq(bmap, &progress);
1468                 bmap = isl_basic_map_gauss(bmap, &progress);
1469                 bmap = remove_duplicate_divs(bmap, &progress);
1470                 bmap = remove_duplicate_constraints(bmap, &progress);
1471         }
1472         return bmap;
1473 }
1474
1475 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1476 {
1477         return (struct isl_basic_set *)
1478                 isl_basic_map_simplify((struct isl_basic_map *)bset);
1479 }
1480
1481 static void dump_term(struct isl_basic_map *bmap,
1482                         isl_int c, int pos, FILE *out)
1483 {
1484         unsigned in = bmap->n_in;
1485         unsigned dim = bmap->n_in + bmap->n_out;
1486         if (!pos)
1487                 isl_int_print(out, c, 0);
1488         else {
1489                 if (!isl_int_is_one(c))
1490                         isl_int_print(out, c, 0);
1491                 if (pos < 1 + bmap->nparam)
1492                         fprintf(out, "p%d", pos - 1);
1493                 else if (pos < 1 + bmap->nparam + in)
1494                         fprintf(out, "i%d", pos - 1 - bmap->nparam);
1495                 else if (pos < 1 + bmap->nparam + dim)
1496                         fprintf(out, "o%d", pos - 1 - bmap->nparam - in);
1497                 else
1498                         fprintf(out, "e%d", pos - 1 - bmap->nparam - dim);
1499         }
1500 }
1501
1502 static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
1503                                 int sign, FILE *out)
1504 {
1505         int i;
1506         int first;
1507         unsigned len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1508         isl_int v;
1509
1510         isl_int_init(v);
1511         for (i = 0, first = 1; i < len; ++i) {
1512                 if (isl_int_sgn(c[i]) * sign <= 0)
1513                         continue;
1514                 if (!first)
1515                         fprintf(out, " + ");
1516                 first = 0;
1517                 isl_int_abs(v, c[i]);
1518                 dump_term(bmap, v, i, out);
1519         }
1520         isl_int_clear(v);
1521         if (first)
1522                 fprintf(out, "0");
1523 }
1524
1525 static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
1526                                 const char *op, FILE *out, int indent)
1527 {
1528         int i;
1529
1530         fprintf(out, "%*s", indent, "");
1531
1532         dump_constraint_sign(bmap, c, 1, out);
1533         fprintf(out, " %s ", op);
1534         dump_constraint_sign(bmap, c, -1, out);
1535
1536         fprintf(out, "\n");
1537
1538         for (i = bmap->n_div; i < bmap->extra; ++i) {
1539                 if (isl_int_is_zero(c[1+bmap->nparam+bmap->n_in+bmap->n_out+i]))
1540                         continue;
1541                 fprintf(out, "%*s", indent, "");
1542                 fprintf(out, "ERROR: unused div coefficient not zero\n");
1543         }
1544 }
1545
1546 static void dump_constraints(struct isl_basic_map *bmap,
1547                                 isl_int **c, unsigned n,
1548                                 const char *op, FILE *out, int indent)
1549 {
1550         int i;
1551
1552         for (i = 0; i < n; ++i)
1553                 dump_constraint(bmap, c[i], op, out, indent);
1554 }
1555
1556 static void dump_affine(struct isl_basic_map *bmap, isl_int *exp, FILE *out)
1557 {
1558         int j;
1559         int first = 1;
1560         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1561
1562         for (j = 0; j < 1 + total; ++j) {
1563                 if (isl_int_is_zero(exp[j]))
1564                         continue;
1565                 if (!first && isl_int_is_pos(exp[j]))
1566                         fprintf(out, "+");
1567                 dump_term(bmap, exp[j], j, out);
1568                 first = 0;
1569         }
1570 }
1571
1572 static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
1573 {
1574         int i;
1575
1576         dump_constraints(bmap, bmap->eq, bmap->n_eq, "=", out, indent);
1577         dump_constraints(bmap, bmap->ineq, bmap->n_ineq, ">=", out, indent);
1578
1579         for (i = 0; i < bmap->n_div; ++i) {
1580                 fprintf(out, "%*s", indent, "");
1581                 fprintf(out, "e%d = [(", i);
1582                 dump_affine(bmap, bmap->div[i]+1, out);
1583                 fprintf(out, ")/");
1584                 isl_int_print(out, bmap->div[i][0], 0);
1585                 fprintf(out, "]\n");
1586         }
1587 }
1588
1589 void isl_basic_set_dump(struct isl_basic_set *bset, FILE *out, int indent)
1590 {
1591         if (!bset) {
1592                 fprintf(out, "null basic set\n");
1593                 return;
1594         }
1595
1596         fprintf(out, "%*s", indent, "");
1597         fprintf(out, "ref: %d, nparam: %d, dim: %d, extra: %d, flags: %x\n",
1598                         bset->ref, bset->nparam, bset->dim, bset->extra,
1599                         bset->flags);
1600         dump((struct isl_basic_map *)bset, out, indent);
1601 }
1602
1603 void isl_basic_map_dump(struct isl_basic_map *bmap, FILE *out, int indent)
1604 {
1605         if (!bmap) {
1606                 fprintf(out, "null basic map\n");
1607                 return;
1608         }
1609
1610         fprintf(out, "%*s", indent, "");
1611         fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d, flags: %x\n",
1612                 bmap->ref,
1613                 bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, bmap->flags);
1614         dump(bmap, out, indent);
1615 }
1616
1617 int isl_inequality_negate(struct isl_basic_map *bmap, unsigned pos)
1618 {
1619         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1620         if (!bmap)
1621                 return -1;
1622         isl_assert(bmap->ctx, pos < bmap->n_ineq, return -1);
1623         isl_seq_neg(bmap->ineq[pos], bmap->ineq[pos], 1 + total);
1624         isl_int_sub_ui(bmap->ineq[pos][0], bmap->ineq[pos][0], 1);
1625         F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1626         return 0;
1627 }
1628
1629 struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
1630                 unsigned nparam, unsigned dim, int n, unsigned flags)
1631 {
1632         struct isl_set *set;
1633
1634         isl_assert(ctx, n >= 0, return NULL);
1635         set = isl_alloc(ctx, struct isl_set,
1636                         sizeof(struct isl_set) +
1637                         n * sizeof(struct isl_basic_set *));
1638         if (!set)
1639                 return NULL;
1640
1641         set->ctx = ctx;
1642         isl_ctx_ref(ctx);
1643         set->ref = 1;
1644         set->size = n;
1645         set->n = 0;
1646         set->nparam = nparam;
1647         set->zero = 0;
1648         set->dim = dim;
1649         set->flags = flags;
1650         return set;
1651 }
1652
1653 struct isl_set *isl_set_dup(struct isl_set *set)
1654 {
1655         int i;
1656         struct isl_set *dup;
1657
1658         if (!set)
1659                 return NULL;
1660
1661         dup = isl_set_alloc(set->ctx, set->nparam, set->dim, set->n, set->flags);
1662         if (!dup)
1663                 return NULL;
1664         for (i = 0; i < set->n; ++i)
1665                 dup = isl_set_add(dup,
1666                                   isl_basic_set_dup(set->p[i]));
1667         return dup;
1668 }
1669
1670 struct isl_set *isl_set_from_basic_set(struct isl_basic_set *bset)
1671 {
1672         struct isl_set *set;
1673
1674         if (!bset)
1675                 return NULL;
1676
1677         set = isl_set_alloc(bset->ctx, bset->nparam, bset->dim, 1, ISL_MAP_DISJOINT);
1678         if (!set) {
1679                 isl_basic_set_free(bset);
1680                 return NULL;
1681         }
1682         return isl_set_add(set, bset);
1683 }
1684
1685 struct isl_map *isl_map_from_basic_map(struct isl_basic_map *bmap)
1686 {
1687         struct isl_map *map;
1688
1689         if (!bmap)
1690                 return NULL;
1691
1692         map = isl_map_alloc(bmap->ctx, bmap->nparam, bmap->n_in, bmap->n_out, 1,
1693                                 ISL_MAP_DISJOINT);
1694         if (!map) {
1695                 isl_basic_map_free(bmap);
1696                 return NULL;
1697         }
1698         return isl_map_add(map, bmap);
1699 }
1700
1701 struct isl_set *isl_set_add(struct isl_set *set, struct isl_basic_set *bset)
1702 {
1703         if (!bset || !set)
1704                 goto error;
1705         isl_assert(set->ctx, set->nparam == bset->nparam, goto error);
1706         isl_assert(set->ctx, set->dim == bset->dim, goto error);
1707         isl_assert(set->ctx, set->n < set->size, goto error);
1708         set->p[set->n] = bset;
1709         set->n++;
1710         return set;
1711 error:
1712         if (set)
1713                 isl_set_free(set);
1714         if (bset)
1715                 isl_basic_set_free(bset);
1716         return NULL;
1717 }
1718
1719 void isl_set_free(struct isl_set *set)
1720 {
1721         int i;
1722
1723         if (!set)
1724                 return;
1725
1726         if (--set->ref > 0)
1727                 return;
1728
1729         isl_ctx_deref(set->ctx);
1730         for (i = 0; i < set->n; ++i)
1731                 isl_basic_set_free(set->p[i]);
1732         free(set);
1733 }
1734
1735 void isl_set_dump(struct isl_set *set, FILE *out, int indent)
1736 {
1737         int i;
1738
1739         if (!set) {
1740                 fprintf(out, "null set\n");
1741                 return;
1742         }
1743
1744         fprintf(out, "%*s", indent, "");
1745         fprintf(out, "ref: %d, n: %d, nparam: %d, dim: %d, flags: %x\n",
1746                         set->ref, set->n, set->nparam, set->dim, set->flags);
1747         for (i = 0; i < set->n; ++i) {
1748                 fprintf(out, "%*s", indent, "");
1749                 fprintf(out, "basic set %d:\n", i);
1750                 isl_basic_set_dump(set->p[i], out, indent+4);
1751         }
1752 }
1753
1754 void isl_map_dump(struct isl_map *map, FILE *out, int indent)
1755 {
1756         int i;
1757
1758         if (!map) {
1759                 fprintf(out, "null map\n");
1760                 return;
1761         }
1762
1763         fprintf(out, "%*s", indent, "");
1764         fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d, flags: %x\n",
1765                         map->ref, map->n, map->nparam, map->n_in, map->n_out,
1766                         map->flags);
1767         for (i = 0; i < map->n; ++i) {
1768                 fprintf(out, "%*s", indent, "");
1769                 fprintf(out, "basic map %d:\n", i);
1770                 isl_basic_map_dump(map->p[i], out, indent+4);
1771         }
1772 }
1773
1774 struct isl_basic_map *isl_basic_map_intersect_domain(
1775                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1776 {
1777         struct isl_basic_map *bmap_domain;
1778
1779         if (!bmap || !bset)
1780                 goto error;
1781
1782         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
1783         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1784
1785         bmap = isl_basic_map_extend(bmap,
1786                         bset->nparam, bmap->n_in, bmap->n_out,
1787                         bset->n_div, bset->n_eq, bset->n_ineq);
1788         if (!bmap)
1789                 goto error;
1790         bmap_domain = isl_basic_map_from_basic_set(bset,
1791                                                 bset->dim, 0);
1792         bmap = add_constraints(bmap, bmap_domain, 0, 0);
1793
1794         bmap = isl_basic_map_simplify(bmap);
1795         return isl_basic_map_finalize(bmap);
1796 error:
1797         isl_basic_map_free(bmap);
1798         isl_basic_set_free(bset);
1799         return NULL;
1800 }
1801
1802 struct isl_basic_map *isl_basic_map_intersect_range(
1803                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1804 {
1805         struct isl_basic_map *bmap_range;
1806
1807         if (!bmap || !bset)
1808                 goto error;
1809
1810         isl_assert(bset->ctx, bset->dim == bmap->n_out, goto error);
1811         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1812
1813         bmap = isl_basic_map_extend(bmap,
1814                         bset->nparam, bmap->n_in, bmap->n_out,
1815                         bset->n_div, bset->n_eq, bset->n_ineq);
1816         if (!bmap)
1817                 goto error;
1818         bmap_range = isl_basic_map_from_basic_set(bset,
1819                                                 0, bset->dim);
1820         bmap = add_constraints(bmap, bmap_range, 0, 0);
1821
1822         bmap = isl_basic_map_simplify(bmap);
1823         return isl_basic_map_finalize(bmap);
1824 error:
1825         isl_basic_map_free(bmap);
1826         isl_basic_set_free(bset);
1827         return NULL;
1828 }
1829
1830 struct isl_basic_map *isl_basic_map_intersect(
1831                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1832 {
1833         if (!bmap1 || !bmap2)
1834                 goto error;
1835
1836         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
1837         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
1838         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
1839
1840         bmap1 = isl_basic_map_extend(bmap1,
1841                         bmap1->nparam, bmap1->n_in, bmap1->n_out,
1842                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
1843         if (!bmap1)
1844                 goto error;
1845         bmap1 = add_constraints(bmap1, bmap2, 0, 0);
1846
1847         bmap1 = isl_basic_map_simplify(bmap1);
1848         return isl_basic_map_finalize(bmap1);
1849 error:
1850         isl_basic_map_free(bmap1);
1851         isl_basic_map_free(bmap2);
1852         return NULL;
1853 }
1854
1855 struct isl_basic_set *isl_basic_set_intersect(
1856                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1857 {
1858         return (struct isl_basic_set *)
1859                 isl_basic_map_intersect(
1860                         (struct isl_basic_map *)bset1,
1861                         (struct isl_basic_map *)bset2);
1862 }
1863
1864 struct isl_map *isl_map_intersect(struct isl_map *map1, struct isl_map *map2)
1865 {
1866         unsigned flags = 0;
1867         struct isl_map *result;
1868         int i, j;
1869
1870         if (!map1 || !map2)
1871                 goto error;
1872
1873         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
1874             F_ISSET(map2, ISL_MAP_DISJOINT))
1875                 FL_SET(flags, ISL_MAP_DISJOINT);
1876
1877         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
1878                                 map1->n * map2->n, flags);
1879         if (!result)
1880                 goto error;
1881         for (i = 0; i < map1->n; ++i)
1882                 for (j = 0; j < map2->n; ++j) {
1883                         struct isl_basic_map *part;
1884                         part = isl_basic_map_intersect(
1885                                     isl_basic_map_copy(map1->p[i]),
1886                                     isl_basic_map_copy(map2->p[j]));
1887                         if (isl_basic_map_is_empty(part))
1888                                 isl_basic_map_free(part);
1889                         else
1890                                 result = isl_map_add(result, part);
1891                         if (!result)
1892                                 goto error;
1893                 }
1894         isl_map_free(map1);
1895         isl_map_free(map2);
1896         return result;
1897 error:
1898         isl_map_free(map1);
1899         isl_map_free(map2);
1900         return NULL;
1901 }
1902
1903 struct isl_set *isl_set_intersect(struct isl_set *set1, struct isl_set *set2)
1904 {
1905         return (struct isl_set *)
1906                 isl_map_intersect((struct isl_map *)set1,
1907                                   (struct isl_map *)set2);
1908 }
1909
1910 struct isl_basic_map *isl_basic_map_reverse(struct isl_basic_map *bmap)
1911 {
1912         struct isl_basic_set *bset;
1913         unsigned in;
1914
1915         if (!bmap)
1916                 return NULL;
1917         bmap = isl_basic_map_cow(bmap);
1918         if (!bmap)
1919                 return NULL;
1920         in = bmap->n_in;
1921         bset = isl_basic_set_from_basic_map(bmap);
1922         bset = isl_basic_set_swap_vars(bset, in);
1923         return isl_basic_map_from_basic_set(bset, bset->dim-in, in);
1924 }
1925
1926 /* Turn final n dimensions into existentially quantified variables.
1927  */
1928 struct isl_basic_set *isl_basic_set_project_out(
1929                 struct isl_basic_set *bset, unsigned n, unsigned flags)
1930 {
1931         int i;
1932         size_t row_size;
1933         isl_int **new_div;
1934         isl_int *old;
1935
1936         if (!bset)
1937                 return NULL;
1938
1939         isl_assert(bset->ctx, n <= bset->dim, goto error);
1940
1941         if (n == 0)
1942                 return bset;
1943
1944         bset = isl_basic_set_cow(bset);
1945
1946         row_size = 1 + bset->nparam + bset->dim + bset->extra;
1947         old = bset->block2.data;
1948         bset->block2 = isl_blk_extend(bset->ctx, bset->block2,
1949                                         (bset->extra + n) * (1 + row_size));
1950         if (!bset->block2.data)
1951                 goto error;
1952         new_div = isl_alloc_array(ctx, isl_int *, bset->extra + n);
1953         if (!new_div)
1954                 goto error;
1955         for (i = 0; i < n; ++i) {
1956                 new_div[i] = bset->block2.data +
1957                                 (bset->extra + i) * (1 + row_size);
1958                 isl_seq_clr(new_div[i], 1 + row_size);
1959         }
1960         for (i = 0; i < bset->extra; ++i)
1961                 new_div[n + i] = bset->block2.data + (bset->div[i] - old);
1962         free(bset->div);
1963         bset->div = new_div;
1964         bset->n_div += n;
1965         bset->extra += n;
1966         bset->dim -= n;
1967         bset = isl_basic_set_simplify(bset);
1968         return isl_basic_set_finalize(bset);
1969 error:
1970         isl_basic_set_free(bset);
1971         return NULL;
1972 }
1973
1974 struct isl_basic_map *isl_basic_map_apply_range(
1975                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1976 {
1977         struct isl_basic_set *bset;
1978         unsigned n_in, n_out;
1979
1980         if (!bmap1 || !bmap2)
1981                 goto error;
1982
1983         isl_assert(bmap->ctx, bmap1->n_out == bmap2->n_in, goto error);
1984         isl_assert(bmap->ctx, bmap1->nparam == bmap2->nparam, goto error);
1985
1986         n_in = bmap1->n_in;
1987         n_out = bmap2->n_out;
1988
1989         bmap2 = isl_basic_map_reverse(bmap2);
1990         if (!bmap2)
1991                 goto error;
1992         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
1993                         bmap1->n_in + bmap2->n_in, bmap1->n_out,
1994                         bmap2->extra, bmap2->n_eq, bmap2->n_ineq);
1995         if (!bmap1)
1996                 goto error;
1997         bmap1 = add_constraints(bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
1998         bmap1 = isl_basic_map_simplify(bmap1);
1999         bset = isl_basic_set_from_basic_map(bmap1);
2000         bset = isl_basic_set_project_out(bset,
2001                                                 bset->dim - (n_in + n_out), 0);
2002         return isl_basic_map_from_basic_set(bset, n_in, n_out);
2003 error:
2004         isl_basic_map_free(bmap1);
2005         isl_basic_map_free(bmap2);
2006         return NULL;
2007 }
2008
2009 struct isl_basic_set *isl_basic_set_apply(
2010                 struct isl_basic_set *bset, struct isl_basic_map *bmap)
2011 {
2012         if (!bset || !bmap)
2013                 goto error;
2014
2015         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
2016
2017         return (struct isl_basic_set *)
2018                 isl_basic_map_apply_range((struct isl_basic_map *)bset, bmap);
2019 error:
2020         isl_basic_set_free(bset);
2021         isl_basic_map_free(bmap);
2022         return NULL;
2023 }
2024
2025 struct isl_basic_map *isl_basic_map_apply_domain(
2026                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
2027 {
2028         if (!bmap1 || !bmap2)
2029                 goto error;
2030
2031         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
2032         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
2033
2034         bmap1 = isl_basic_map_reverse(bmap1);
2035         bmap1 = isl_basic_map_apply_range(bmap1, bmap2);
2036         return isl_basic_map_reverse(bmap1);
2037 error:
2038         isl_basic_map_free(bmap1);
2039         isl_basic_map_free(bmap2);
2040         return NULL;
2041 }
2042
2043 static struct isl_basic_map *var_equal(struct isl_ctx *ctx,
2044                 struct isl_basic_map *bmap, unsigned pos)
2045 {
2046         int i;
2047         i = isl_basic_map_alloc_equality(bmap);
2048         if (i < 0)
2049                 goto error;
2050         isl_seq_clr(bmap->eq[i],
2051                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2052         isl_int_set_si(bmap->eq[i][1+bmap->nparam+pos], -1);
2053         isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+pos], 1);
2054         return isl_basic_map_finalize(bmap);
2055 error:
2056         isl_basic_map_free(bmap);
2057         return NULL;
2058 }
2059
2060 static struct isl_basic_map *var_more(struct isl_ctx *ctx,
2061                 struct isl_basic_map *bmap, unsigned pos)
2062 {
2063         int i;
2064         i = isl_basic_map_alloc_inequality(bmap);
2065         if (i < 0)
2066                 goto error;
2067         isl_seq_clr(bmap->ineq[i],
2068                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2069         isl_int_set_si(bmap->ineq[i][0], -1);
2070         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], -1);
2071         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], 1);
2072         return isl_basic_map_finalize(bmap);
2073 error:
2074         isl_basic_map_free(bmap);
2075         return NULL;
2076 }
2077
2078 static struct isl_basic_map *var_less(struct isl_ctx *ctx,
2079                 struct isl_basic_map *bmap, unsigned pos)
2080 {
2081         int i;
2082         i = isl_basic_map_alloc_inequality(bmap);
2083         if (i < 0)
2084                 goto error;
2085         isl_seq_clr(bmap->ineq[i],
2086                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2087         isl_int_set_si(bmap->ineq[i][0], -1);
2088         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], 1);
2089         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], -1);
2090         return isl_basic_map_finalize(bmap);
2091 error:
2092         isl_basic_map_free(bmap);
2093         return NULL;
2094 }
2095
2096 struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
2097                 unsigned nparam, unsigned in, unsigned out, unsigned n_equal)
2098 {
2099         int i;
2100         struct isl_basic_map *bmap;
2101         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, n_equal, 0);
2102         if (!bmap)
2103                 return NULL;
2104         for (i = 0; i < n_equal && bmap; ++i)
2105                 bmap = var_equal(ctx, bmap, i);
2106         return isl_basic_map_finalize(bmap);
2107 }
2108
2109 struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
2110                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
2111 {
2112         int i;
2113         struct isl_basic_map *bmap;
2114         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
2115         if (!bmap)
2116                 return NULL;
2117         for (i = 0; i < pos && bmap; ++i)
2118                 bmap = var_equal(ctx, bmap, i);
2119         if (bmap)
2120                 bmap = var_less(ctx, bmap, pos);
2121         return isl_basic_map_finalize(bmap);
2122 }
2123
2124 struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
2125                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
2126 {
2127         int i;
2128         struct isl_basic_map *bmap;
2129         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
2130         if (!bmap)
2131                 return NULL;
2132         for (i = 0; i < pos && bmap; ++i)
2133                 bmap = var_equal(ctx, bmap, i);
2134         if (bmap)
2135                 bmap = var_more(ctx, bmap, pos);
2136         return isl_basic_map_finalize(bmap);
2137 }
2138
2139 struct isl_basic_map *isl_basic_map_from_basic_set(
2140                 struct isl_basic_set *bset, unsigned n_in, unsigned n_out)
2141 {
2142         struct isl_basic_map *bmap;
2143
2144         bset = isl_basic_set_cow(bset);
2145         if (!bset)
2146                 return NULL;
2147
2148         isl_assert(ctx, bset->dim == n_in + n_out, goto error);
2149         bmap = (struct isl_basic_map *) bset;
2150         bmap->n_in = n_in;
2151         bmap->n_out = n_out;
2152         return isl_basic_map_finalize(bmap);
2153 error:
2154         isl_basic_set_free(bset);
2155         return NULL;
2156 }
2157
2158 struct isl_basic_set *isl_basic_set_from_basic_map(struct isl_basic_map *bmap)
2159 {
2160         if (!bmap)
2161                 goto error;
2162         if (bmap->n_in == 0)
2163                 return (struct isl_basic_set *)bmap;
2164         bmap = isl_basic_map_cow(bmap);
2165         if (!bmap)
2166                 goto error;
2167         bmap->n_out += bmap->n_in;
2168         bmap->n_in = 0;
2169         bmap = isl_basic_map_finalize(bmap);
2170         return (struct isl_basic_set *)bmap;
2171 error:
2172         return NULL;
2173 }
2174
2175 /* For a div d = floor(f/m), add the constraints
2176  *
2177  *              f - m d >= 0
2178  *              -(f-(n-1)) + m d >= 0
2179  *
2180  * Note that the second constraint is the negation of
2181  *
2182  *              f - m d >= n
2183  */
2184 static int add_div_constraints(struct isl_basic_map *bmap, unsigned div)
2185 {
2186         int i, j;
2187         unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2188         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
2189
2190         i = isl_basic_map_alloc_inequality(bmap);
2191         if (i < 0)
2192                 return -1;
2193         isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
2194         isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
2195
2196         j = isl_basic_map_alloc_inequality(bmap);
2197         if (j < 0)
2198                 return -1;
2199         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
2200         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
2201         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
2202         return j;
2203 }
2204
2205 struct isl_basic_set *isl_basic_map_underlying_set(
2206                 struct isl_basic_map *bmap)
2207 {
2208         if (!bmap)
2209                 goto error;
2210         if (bmap->nparam == 0 && bmap->n_in == 0 && bmap->n_div == 0)
2211                 return (struct isl_basic_set *)bmap;
2212         bmap = isl_basic_map_cow(bmap);
2213         if (!bmap)
2214                 goto error;
2215         bmap->n_out += bmap->nparam + bmap->n_in + bmap->n_div;
2216         bmap->nparam = 0;
2217         bmap->n_in = 0;
2218         bmap->extra -= bmap->n_div;
2219         bmap->n_div = 0;
2220         bmap = isl_basic_map_finalize(bmap);
2221         return (struct isl_basic_set *)bmap;
2222 error:
2223         return NULL;
2224 }
2225
2226 struct isl_basic_map *isl_basic_map_overlying_set(
2227         struct isl_basic_set *bset, struct isl_basic_map *like)
2228 {
2229         struct isl_basic_map *bmap;
2230         struct isl_ctx *ctx;
2231         unsigned total;
2232         int i, k;
2233
2234         if (!bset || !like)
2235                 goto error;
2236         ctx = bset->ctx;
2237         isl_assert(ctx, bset->dim ==
2238                         like->nparam + like->n_in + like->n_out + like->n_div,
2239                         goto error);
2240         if (like->nparam == 0 && like->n_in == 0 && like->n_div == 0) {
2241                 isl_basic_map_free(like);
2242                 return (struct isl_basic_map *)bset;
2243         }
2244         bset = isl_basic_set_cow(bset);
2245         if (!bset)
2246                 goto error;
2247         total = bset->dim + bset->extra;
2248         bmap = (struct isl_basic_map *)bset;
2249         bmap->nparam = like->nparam;
2250         bmap->n_in = like->n_in;
2251         bmap->n_out = like->n_out;
2252         bmap->extra += like->n_div;
2253         if (bmap->extra) {
2254                 unsigned ltotal;
2255                 ltotal = total - bmap->extra + like->extra;
2256                 if (ltotal > total)
2257                         ltotal = total;
2258                 bmap->block2 = isl_blk_extend(ctx, bmap->block2,
2259                                         bmap->extra * (1 + 1 + total));
2260                 if (isl_blk_is_error(bmap->block2))
2261                         goto error;
2262                 bmap->div = isl_realloc_array(ctx, bmap->div, isl_int *,
2263                                                 bmap->extra);
2264                 if (!bmap->div)
2265                         goto error;
2266                 bmap = isl_basic_map_extend(bmap, bmap->nparam,
2267                         bmap->n_in, bmap->n_out, 0, 0, 2 * like->n_div);
2268                 for (i = 0; i < like->n_div; ++i) {
2269                         k = isl_basic_map_alloc_div(bmap);
2270                         if (k < 0)
2271                                 goto error;
2272                         isl_seq_cpy(bmap->div[k], like->div[i], 1 + 1 + ltotal);
2273                         isl_seq_clr(bmap->div[k]+1+1+ltotal, total - ltotal);
2274                         if (add_div_constraints(bmap, k) < 0)
2275                                 goto error;
2276                 }
2277         }
2278         isl_basic_map_free(like);
2279         bmap = isl_basic_map_finalize(bmap);
2280         return bmap;
2281 error:
2282         isl_basic_map_free(like);
2283         isl_basic_set_free(bset);
2284         return NULL;
2285 }
2286
2287 struct isl_basic_set *isl_basic_set_from_underlying_set(
2288         struct isl_basic_set *bset, struct isl_basic_set *like)
2289 {
2290         return (struct isl_basic_set *)
2291                 isl_basic_map_overlying_set(bset, (struct isl_basic_map *)like);
2292 }
2293
2294 struct isl_set *isl_set_from_underlying_set(
2295         struct isl_set *set, struct isl_basic_set *like)
2296 {
2297         int i;
2298
2299         if (!set || !like)
2300                 goto error;
2301         isl_assert(set->ctx, set->dim == like->nparam + like->dim + like->n_div,
2302                     goto error);
2303         if (like->nparam == 0 && like->n_div == 0) {
2304                 isl_basic_set_free(like);
2305                 return set;
2306         }
2307         set = isl_set_cow(set);
2308         if (!set)
2309                 goto error;
2310         for (i = 0; i < set->n; ++i) {
2311                 set->p[i] = isl_basic_set_from_underlying_set(set->p[i],
2312                                                       isl_basic_set_copy(like));
2313                 if (!set->p[i])
2314                         goto error;
2315         }
2316         set->nparam = like->nparam;
2317         set->dim = like->dim;
2318         isl_basic_set_free(like);
2319         return set;
2320 error:
2321         isl_basic_set_free(like);
2322         isl_set_free(set);
2323         return NULL;
2324 }
2325
2326 struct isl_set *isl_map_underlying_set(struct isl_map *map)
2327 {
2328         int i;
2329
2330         map = isl_map_align_divs(map);
2331         map = isl_map_cow(map);
2332         if (!map)
2333                 return NULL;
2334
2335         for (i = 0; i < map->n; ++i) {
2336                 map->p[i] = (struct isl_basic_map *)
2337                                 isl_basic_map_underlying_set(map->p[i]);
2338                 if (!map->p[i])
2339                         goto error;
2340         }
2341         if (map->n == 0)
2342                 map->n_out += map->nparam + map->n_in;
2343         else
2344                 map->n_out = map->p[0]->n_out;
2345         map->nparam = 0;
2346         map->n_in = 0;
2347         return (struct isl_set *)map;
2348 error:
2349         isl_map_free(map);
2350         return NULL;
2351 }
2352
2353 struct isl_set *isl_set_to_underlying_set(struct isl_set *set)
2354 {
2355         return (struct isl_set *)isl_map_underlying_set((struct isl_map *)set);
2356 }
2357
2358 struct isl_basic_set *isl_basic_map_domain(struct isl_basic_map *bmap)
2359 {
2360         struct isl_basic_set *domain;
2361         unsigned n_out;
2362         if (!bmap)
2363                 return NULL;
2364         n_out = bmap->n_out;
2365         domain = isl_basic_set_from_basic_map(bmap);
2366         return isl_basic_set_project_out(domain, n_out, 0);
2367 }
2368
2369 struct isl_basic_set *isl_basic_map_range(struct isl_basic_map *bmap)
2370 {
2371         return isl_basic_map_domain(isl_basic_map_reverse(bmap));
2372 }
2373
2374 struct isl_set *isl_map_range(struct isl_map *map)
2375 {
2376         int i;
2377         struct isl_set *set;
2378
2379         if (!map)
2380                 goto error;
2381         map = isl_map_cow(map);
2382         if (!map)
2383                 goto error;
2384
2385         set = (struct isl_set *) map;
2386         set->zero = 0;
2387         for (i = 0; i < map->n; ++i) {
2388                 set->p[i] = isl_basic_map_range(map->p[i]);
2389                 if (!set->p[i])
2390                         goto error;
2391         }
2392         F_CLR(set, ISL_MAP_DISJOINT);
2393         F_CLR(set, ISL_SET_NORMALIZED);
2394         return set;
2395 error:
2396         isl_map_free(map);
2397         return NULL;
2398 }
2399
2400 struct isl_map *isl_map_from_set(struct isl_set *set,
2401                 unsigned n_in, unsigned n_out)
2402 {
2403         int i;
2404         struct isl_map *map = NULL;
2405
2406         if (!set)
2407                 return NULL;
2408         isl_assert(set->ctx, set->dim == n_in + n_out, goto error);
2409         set = isl_set_cow(set);
2410         if (!set)
2411                 return NULL;
2412         map = (struct isl_map *)set;
2413         for (i = 0; i < set->n; ++i) {
2414                 map->p[i] = isl_basic_map_from_basic_set(
2415                                 set->p[i], n_in, n_out);
2416                 if (!map->p[i])
2417                         goto error;
2418         }
2419         map->n_in = n_in;
2420         map->n_out = n_out;
2421         return map;
2422 error:
2423         isl_set_free(set);
2424         return NULL;
2425 }
2426
2427 struct isl_set *isl_set_from_map(struct isl_map *map)
2428 {
2429         int i;
2430         struct isl_set *set = NULL;
2431
2432         if (!map)
2433                 return NULL;
2434         map = isl_map_cow(map);
2435         if (!map)
2436                 return NULL;
2437         map->n_out += map->n_in;
2438         map->n_in = 0;
2439         set = (struct isl_set *)map;
2440         for (i = 0; i < map->n; ++i) {
2441                 set->p[i] = isl_basic_set_from_basic_map(map->p[i]);
2442                 if (!set->p[i])
2443                         goto error;
2444         }
2445         return set;
2446 error:
2447         isl_map_free(map);
2448         return NULL;
2449 }
2450
2451 struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
2452                 unsigned nparam, unsigned in, unsigned out, int n,
2453                 unsigned flags)
2454 {
2455         struct isl_map *map;
2456
2457         map = isl_alloc(ctx, struct isl_map,
2458                         sizeof(struct isl_map) +
2459                         n * sizeof(struct isl_basic_map *));
2460         if (!map)
2461                 return NULL;
2462
2463         map->ctx = ctx;
2464         isl_ctx_ref(ctx);
2465         map->ref = 1;
2466         map->size = n;
2467         map->n = 0;
2468         map->nparam = nparam;
2469         map->n_in = in;
2470         map->n_out = out;
2471         map->flags = flags;
2472         return map;
2473 }
2474
2475 struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
2476                 unsigned nparam, unsigned in, unsigned out)
2477 {
2478         struct isl_basic_map *bmap;
2479         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 1, 0);
2480         bmap = isl_basic_map_set_to_empty(bmap);
2481         return bmap;
2482 }
2483
2484 struct isl_basic_set *isl_basic_set_empty(struct isl_ctx *ctx,
2485                 unsigned nparam, unsigned dim)
2486 {
2487         struct isl_basic_set *bset;
2488         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 1, 0);
2489         bset = isl_basic_set_set_to_empty(bset);
2490         return bset;
2491 }
2492
2493 struct isl_basic_map *isl_basic_map_universe(struct isl_ctx *ctx,
2494                 unsigned nparam, unsigned in, unsigned out)
2495 {
2496         struct isl_basic_map *bmap;
2497         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 0, 0);
2498         return bmap;
2499 }
2500
2501 struct isl_basic_set *isl_basic_set_universe(struct isl_ctx *ctx,
2502                 unsigned nparam, unsigned dim)
2503 {
2504         struct isl_basic_set *bset;
2505         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 0, 0);
2506         return bset;
2507 }
2508
2509 struct isl_map *isl_map_empty(struct isl_ctx *ctx,
2510                 unsigned nparam, unsigned in, unsigned out)
2511 {
2512         return isl_map_alloc(ctx, nparam, in, out, 0, ISL_MAP_DISJOINT);
2513 }
2514
2515 struct isl_set *isl_set_empty(struct isl_ctx *ctx,
2516                 unsigned nparam, unsigned dim)
2517 {
2518         return isl_set_alloc(ctx, nparam, dim, 0, ISL_MAP_DISJOINT);
2519 }
2520
2521 struct isl_map *isl_map_dup(struct isl_map *map)
2522 {
2523         int i;
2524         struct isl_map *dup;
2525
2526         if (!map)
2527                 return NULL;
2528         dup = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out, map->n,
2529                                 map->flags);
2530         for (i = 0; i < map->n; ++i)
2531                 dup = isl_map_add(dup,
2532                                   isl_basic_map_dup(map->p[i]));
2533         return dup;
2534 }
2535
2536 struct isl_map *isl_map_add(struct isl_map *map, struct isl_basic_map *bmap)
2537 {
2538         if (!bmap || !map)
2539                 goto error;
2540         isl_assert(map->ctx, map->nparam == bmap->nparam, goto error);
2541         isl_assert(map->ctx, map->n_in == bmap->n_in, goto error);
2542         isl_assert(map->ctx, map->n_out == bmap->n_out, goto error);
2543         isl_assert(map->ctx, map->n < map->size, goto error);
2544         map->p[map->n] = bmap;
2545         map->n++;
2546         F_CLR(map, ISL_MAP_NORMALIZED);
2547         return map;
2548 error:
2549         if (map)
2550                 isl_map_free(map);
2551         if (bmap)
2552                 isl_basic_map_free(bmap);
2553         return NULL;
2554 }
2555
2556 void isl_map_free(struct isl_map *map)
2557 {
2558         int i;
2559
2560         if (!map)
2561                 return;
2562
2563         if (--map->ref > 0)
2564                 return;
2565
2566         isl_ctx_deref(map->ctx);
2567         for (i = 0; i < map->n; ++i)
2568                 isl_basic_map_free(map->p[i]);
2569         free(map);
2570 }
2571
2572 struct isl_map *isl_map_extend(struct isl_map *base,
2573                 unsigned nparam, unsigned n_in, unsigned n_out)
2574 {
2575         int i;
2576
2577         base = isl_map_cow(base);
2578         if (!base)
2579                 return NULL;
2580
2581         isl_assert(base->ctx, base->nparam <= nparam, goto error);
2582         isl_assert(base->ctx, base->n_in <= n_in, goto error);
2583         isl_assert(base->ctx, base->n_out <= n_out, goto error);
2584         base->nparam = nparam;
2585         base->n_in = n_in;
2586         base->n_out = n_out;
2587         for (i = 0; i < base->n; ++i) {
2588                 base->p[i] = isl_basic_map_extend(base->p[i],
2589                                 nparam, n_in, n_out, 0, 0, 0);
2590                 if (!base->p[i])
2591                         goto error;
2592         }
2593         return base;
2594 error:
2595         isl_map_free(base);
2596         return NULL;
2597 }
2598
2599 struct isl_set *isl_set_extend(struct isl_set *base,
2600                 unsigned nparam, unsigned dim)
2601 {
2602         return (struct isl_set *)isl_map_extend((struct isl_map *)base,
2603                                                         nparam, 0, dim);
2604 }
2605
2606 static struct isl_basic_map *isl_basic_map_fix_var(struct isl_basic_map *bmap,
2607                 unsigned var, int value)
2608 {
2609         int j;
2610
2611         bmap = isl_basic_map_cow(bmap);
2612         bmap = isl_basic_map_extend(bmap,
2613                         bmap->nparam, bmap->n_in, bmap->n_out, 0, 1, 0);
2614         j = isl_basic_map_alloc_equality(bmap);
2615         if (j < 0)
2616                 goto error;
2617         isl_seq_clr(bmap->eq[j],
2618                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2619         isl_int_set_si(bmap->eq[j][1 + var], -1);
2620         isl_int_set_si(bmap->eq[j][0], value);
2621         bmap = isl_basic_map_simplify(bmap);
2622         return isl_basic_map_finalize(bmap);
2623 error:
2624         isl_basic_map_free(bmap);
2625         return NULL;
2626 }
2627
2628 struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_basic_map *bmap,
2629                 unsigned input, int value)
2630 {
2631         if (!bmap)
2632                 return NULL;
2633         isl_assert(bmap->ctx, input < bmap->n_in, goto error);
2634         return isl_basic_map_fix_var(bmap, bmap->nparam + input, value);
2635 error:
2636         isl_basic_map_free(bmap);
2637         return NULL;
2638 }
2639
2640 struct isl_basic_set *isl_basic_set_fix_dim_si(struct isl_basic_set *bset,
2641                 unsigned dim, int value)
2642 {
2643         if (!bset)
2644                 return NULL;
2645         isl_assert(bset->ctx, dim < bset->dim, goto error);
2646         return (struct isl_basic_set *)
2647                 isl_basic_map_fix_var((struct isl_basic_map *)bset,
2648                                                 bset->nparam + dim, value);
2649 error:
2650         isl_basic_set_free(bset);
2651         return NULL;
2652 }
2653
2654 struct isl_map *isl_map_fix_input_si(struct isl_map *map,
2655                 unsigned input, int value)
2656 {
2657         int i;
2658
2659         map = isl_map_cow(map);
2660         if (!map)
2661                 return NULL;
2662
2663         isl_assert(ctx, input < map->n_in, goto error);
2664         for (i = 0; i < map->n; ++i) {
2665                 map->p[i] = isl_basic_map_fix_input_si(map->p[i], input, value);
2666                 if (!map->p[i])
2667                         goto error;
2668         }
2669         F_CLR(map, ISL_MAP_NORMALIZED);
2670         return map;
2671 error:
2672         isl_map_free(map);
2673         return NULL;
2674 }
2675
2676 struct isl_set *isl_set_fix_dim_si(struct isl_set *set, unsigned dim, int value)
2677 {
2678         int i;
2679
2680         set = isl_set_cow(set);
2681         if (!set)
2682                 return NULL;
2683
2684         isl_assert(set->ctx, dim < set->dim, goto error);
2685         for (i = 0; i < set->n; ++i) {
2686                 set->p[i] = isl_basic_set_fix_dim_si(set->p[i], dim, value);
2687                 if (!set->p[i])
2688                         goto error;
2689         }
2690         return set;
2691 error:
2692         isl_set_free(set);
2693         return NULL;
2694 }
2695
2696 struct isl_basic_set *isl_basic_set_lower_bound_dim(struct isl_basic_set *bset,
2697         unsigned dim, isl_int value)
2698 {
2699         int j;
2700
2701         bset = isl_basic_set_cow(bset);
2702         bset = isl_basic_set_extend(bset, bset->nparam, bset->dim, 0, 0, 1);
2703         j = isl_basic_set_alloc_inequality(bset);
2704         if (j < 0)
2705                 goto error;
2706         isl_seq_clr(bset->ineq[j], 1 + bset->nparam + bset->dim + bset->extra);
2707         isl_int_set_si(bset->ineq[j][1 + bset->nparam + dim], 1);
2708         isl_int_neg(bset->ineq[j][0], value);
2709         bset = isl_basic_set_simplify(bset);
2710         return isl_basic_set_finalize(bset);
2711 error:
2712         isl_basic_set_free(bset);
2713         return NULL;
2714 }
2715
2716 struct isl_set *isl_set_lower_bound_dim(struct isl_set *set, unsigned dim,
2717                                         isl_int value)
2718 {
2719         int i;
2720
2721         set = isl_set_cow(set);
2722         if (!set)
2723                 return NULL;
2724
2725         isl_assert(set->ctx, dim < set->dim, goto error);
2726         for (i = 0; i < set->n; ++i) {
2727                 set->p[i] = isl_basic_set_lower_bound_dim(set->p[i], dim, value);
2728                 if (!set->p[i])
2729                         goto error;
2730         }
2731         return set;
2732 error:
2733         isl_set_free(set);
2734         return NULL;
2735 }
2736
2737 struct isl_map *isl_map_reverse(struct isl_map *map)
2738 {
2739         int i;
2740         unsigned t;
2741
2742         map = isl_map_cow(map);
2743         if (!map)
2744                 return NULL;
2745
2746         t = map->n_in;
2747         map->n_in = map->n_out;
2748         map->n_out = t;
2749         for (i = 0; i < map->n; ++i) {
2750                 map->p[i] = isl_basic_map_reverse(map->p[i]);
2751                 if (!map->p[i])
2752                         goto error;
2753         }
2754         F_CLR(map, ISL_MAP_NORMALIZED);
2755         return map;
2756 error:
2757         isl_map_free(map);
2758         return NULL;
2759 }
2760
2761 struct isl_map *isl_basic_map_lexmax(
2762                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2763                 struct isl_set **empty)
2764 {
2765         return isl_pip_basic_map_lexmax(bmap, dom, empty);
2766 }
2767
2768 struct isl_map *isl_basic_map_lexmin(
2769                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2770                 struct isl_set **empty)
2771 {
2772         return isl_pip_basic_map_lexmin(bmap, dom, empty);
2773 }
2774
2775 struct isl_set *isl_basic_set_lexmin(struct isl_basic_set *bset)
2776 {
2777         struct isl_basic_map *bmap = NULL;
2778         struct isl_basic_set *dom = NULL;
2779         struct isl_map *min;
2780
2781         if (!bset)
2782                 goto error;
2783         bmap = isl_basic_map_from_basic_set(bset, 0, bset->dim);
2784         if (!bmap)
2785                 goto error;
2786         dom = isl_basic_set_alloc(bmap->ctx, bmap->nparam, 0, 0, 0, 0);
2787         if (!dom)
2788                 goto error;
2789         min = isl_basic_map_lexmin(bmap, dom, NULL);
2790         return isl_map_range(min);
2791 error:
2792         isl_basic_map_free(bmap);
2793         return NULL;
2794 }
2795
2796 struct isl_map *isl_basic_map_compute_divs(struct isl_basic_map *bmap)
2797 {
2798         if (!bmap)
2799                 return NULL;
2800         if (bmap->n_div == 0)
2801                 return isl_map_from_basic_map(bmap);
2802         return isl_pip_basic_map_compute_divs(bmap);
2803 }
2804
2805 struct isl_map *isl_map_compute_divs(struct isl_map *map)
2806 {
2807         int i;
2808         struct isl_map *res;
2809
2810         if (!map)
2811                 return NULL;
2812         if (map->n == 0)
2813                 return map;
2814         res = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[0]));
2815         for (i = 1 ; i < map->n; ++i) {
2816                 struct isl_map *r2;
2817                 r2 = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[i]));
2818                 if (F_ISSET(map, ISL_MAP_DISJOINT))
2819                         res = isl_map_union_disjoint(res, r2);
2820                 else
2821                         res = isl_map_union(res, r2);
2822         }
2823         isl_map_free(map);
2824
2825         return res;
2826 }
2827
2828 struct isl_set *isl_basic_set_compute_divs(struct isl_basic_set *bset)
2829 {
2830         return (struct isl_set *)
2831                 isl_basic_map_compute_divs((struct isl_basic_map *)bset);
2832 }
2833
2834 struct isl_set *isl_set_compute_divs(struct isl_set *set)
2835 {
2836         return (struct isl_set *)
2837                 isl_map_compute_divs((struct isl_map *)set);
2838 }
2839
2840 struct isl_set *isl_map_domain(struct isl_map *map)
2841 {
2842         int i;
2843         struct isl_set *set;
2844
2845         if (!map)
2846                 goto error;
2847
2848         map = isl_map_cow(map);
2849         if (!map)
2850                 return NULL;
2851
2852         set = (struct isl_set *)map;
2853         set->dim = map->n_in;
2854         set->zero = 0;
2855         for (i = 0; i < map->n; ++i) {
2856                 set->p[i] = isl_basic_map_domain(map->p[i]);
2857                 if (!set->p[i])
2858                         goto error;
2859         }
2860         F_CLR(set, ISL_MAP_DISJOINT);
2861         F_CLR(set, ISL_SET_NORMALIZED);
2862         return set;
2863 error:
2864         isl_map_free(map);
2865         return NULL;
2866 }
2867
2868 struct isl_map *isl_map_union_disjoint(
2869                         struct isl_map *map1, struct isl_map *map2)
2870 {
2871         int i;
2872         unsigned flags = 0;
2873         struct isl_map *map = NULL;
2874
2875         if (!map1 || !map2)
2876                 goto error;
2877
2878         if (map1->n == 0) {
2879                 isl_map_free(map1);
2880                 return map2;
2881         }
2882         if (map2->n == 0) {
2883                 isl_map_free(map2);
2884                 return map1;
2885         }
2886
2887         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
2888         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
2889         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
2890
2891         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
2892             F_ISSET(map2, ISL_MAP_DISJOINT))
2893                 FL_SET(flags, ISL_MAP_DISJOINT);
2894
2895         map = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
2896                                 map1->n + map2->n, flags);
2897         if (!map)
2898                 goto error;
2899         for (i = 0; i < map1->n; ++i) {
2900                 map = isl_map_add(map,
2901                                   isl_basic_map_copy(map1->p[i]));
2902                 if (!map)
2903                         goto error;
2904         }
2905         for (i = 0; i < map2->n; ++i) {
2906                 map = isl_map_add(map,
2907                                   isl_basic_map_copy(map2->p[i]));
2908                 if (!map)
2909                         goto error;
2910         }
2911         isl_map_free(map1);
2912         isl_map_free(map2);
2913         return map;
2914 error:
2915         isl_map_free(map);
2916         isl_map_free(map1);
2917         isl_map_free(map2);
2918         return NULL;
2919 }
2920
2921 struct isl_map *isl_map_union(struct isl_map *map1, struct isl_map *map2)
2922 {
2923         map1 = isl_map_union_disjoint(map1, map2);
2924         if (!map1)
2925                 return NULL;
2926         if (map1->n > 1)
2927                 F_CLR(map1, ISL_MAP_DISJOINT);
2928         return map1;
2929 }
2930
2931 struct isl_set *isl_set_union_disjoint(
2932                         struct isl_set *set1, struct isl_set *set2)
2933 {
2934         return (struct isl_set *)
2935                 isl_map_union_disjoint(
2936                         (struct isl_map *)set1, (struct isl_map *)set2);
2937 }
2938
2939 struct isl_set *isl_set_union(struct isl_set *set1, struct isl_set *set2)
2940 {
2941         return (struct isl_set *)
2942                 isl_map_union((struct isl_map *)set1, (struct isl_map *)set2);
2943 }
2944
2945 struct isl_map *isl_map_intersect_range(
2946                 struct isl_map *map, struct isl_set *set)
2947 {
2948         unsigned flags = 0;
2949         struct isl_map *result;
2950         int i, j;
2951
2952         if (!map || !set)
2953                 goto error;
2954
2955         if (F_ISSET(map, ISL_MAP_DISJOINT) &&
2956             F_ISSET(set, ISL_MAP_DISJOINT))
2957                 FL_SET(flags, ISL_MAP_DISJOINT);
2958
2959         result = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
2960                                 map->n * set->n, flags);
2961         if (!result)
2962                 goto error;
2963         for (i = 0; i < map->n; ++i)
2964                 for (j = 0; j < set->n; ++j) {
2965                         result = isl_map_add(result,
2966                             isl_basic_map_intersect_range(
2967                                 isl_basic_map_copy(map->p[i]),
2968                                 isl_basic_set_copy(set->p[j])));
2969                         if (!result)
2970                                 goto error;
2971                 }
2972         isl_map_free(map);
2973         isl_set_free(set);
2974         return result;
2975 error:
2976         isl_map_free(map);
2977         isl_set_free(set);
2978         return NULL;
2979 }
2980
2981 struct isl_map *isl_map_intersect_domain(
2982                 struct isl_map *map, struct isl_set *set)
2983 {
2984         return isl_map_reverse(
2985                 isl_map_intersect_range(isl_map_reverse(map), set));
2986 }
2987
2988 struct isl_map *isl_map_apply_domain(
2989                 struct isl_map *map1, struct isl_map *map2)
2990 {
2991         if (!map1 || !map2)
2992                 goto error;
2993         map1 = isl_map_reverse(map1);
2994         map1 = isl_map_apply_range(map1, map2);
2995         return isl_map_reverse(map1);
2996 error:
2997         isl_map_free(map1);
2998         isl_map_free(map2);
2999         return NULL;
3000 }
3001
3002 struct isl_map *isl_map_apply_range(
3003                 struct isl_map *map1, struct isl_map *map2)
3004 {
3005         struct isl_map *result;
3006         int i, j;
3007
3008         if (!map1 || !map2)
3009                 goto error;
3010
3011         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
3012         isl_assert(map1->ctx, map1->n_out == map2->n_in, goto error);
3013
3014         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map2->n_out,
3015                                 map1->n * map2->n, 0);
3016         if (!result)
3017                 goto error;
3018         for (i = 0; i < map1->n; ++i)
3019                 for (j = 0; j < map2->n; ++j) {
3020                         result = isl_map_add(result,
3021                             isl_basic_map_apply_range(
3022                                 isl_basic_map_copy(map1->p[i]),
3023                                 isl_basic_map_copy(map2->p[j])));
3024                         if (!result)
3025                                 goto error;
3026                 }
3027         isl_map_free(map1);
3028         isl_map_free(map2);
3029         if (result && result->n <= 1)
3030                 F_SET(result, ISL_MAP_DISJOINT);
3031         return result;
3032 error:
3033         isl_map_free(map1);
3034         isl_map_free(map2);
3035         return NULL;
3036 }
3037
3038 /*
3039  * returns range - domain
3040  */
3041 struct isl_basic_set *isl_basic_map_deltas(struct isl_basic_map *bmap)
3042 {
3043         struct isl_basic_set *bset;
3044         unsigned dim;
3045         int i;
3046
3047         if (!bmap)
3048                 goto error;
3049         isl_assert(bmap->ctx, bmap->n_in == bmap->n_out, goto error);
3050         dim = bmap->n_in;
3051         bset = isl_basic_set_from_basic_map(bmap);
3052         bset = isl_basic_set_extend(bset, bset->nparam, 3*dim, 0,
3053                                         dim, 0);
3054         bset = isl_basic_set_swap_vars(bset, 2*dim);
3055         for (i = 0; i < dim; ++i) {
3056                 int j = isl_basic_map_alloc_equality(
3057                                             (struct isl_basic_map *)bset);
3058                 if (j < 0)
3059                         goto error;
3060                 isl_seq_clr(bset->eq[j],
3061                             1 + bset->nparam + bset->dim + bset->extra);
3062                 isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
3063                 isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
3064                 isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
3065         }
3066         return isl_basic_set_project_out(bset, 2*dim, 0);
3067 error:
3068         isl_basic_map_free(bmap);
3069         return NULL;
3070 }
3071
3072 /*
3073  * returns range - domain
3074  */
3075 struct isl_set *isl_map_deltas(struct isl_map *map)
3076 {
3077         int i;
3078         struct isl_set *result;
3079
3080         if (!map)
3081                 return NULL;
3082
3083         isl_assert(map->ctx, map->n_in == map->n_out, goto error);
3084         result = isl_set_alloc(map->ctx, map->nparam, map->n_in, map->n, map->flags);
3085         if (!result)
3086                 goto error;
3087         for (i = 0; i < map->n; ++i)
3088                 result = isl_set_add(result,
3089                           isl_basic_map_deltas(isl_basic_map_copy(map->p[i])));
3090         isl_map_free(map);
3091         return result;
3092 error:
3093         isl_map_free(map);
3094         return NULL;
3095 }
3096
3097 /* If the only constraints a div d=floor(f/m)
3098  * appears in are its two defining constraints
3099  *
3100  *      f - m d >=0
3101  *      -(f - (m - 1)) + m d >= 0
3102  *
3103  * then it can safely be removed.
3104  */
3105 static int div_is_redundant(struct isl_basic_map *bmap, int div)
3106 {
3107         int i;
3108         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
3109
3110         for (i = 0; i < bmap->n_eq; ++i)
3111                 if (!isl_int_is_zero(bmap->eq[i][pos]))
3112                         return 0;
3113
3114         for (i = 0; i < bmap->n_ineq; ++i) {
3115                 if (isl_int_is_zero(bmap->ineq[i][pos]))
3116                         continue;
3117                 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
3118                         int neg;
3119                         isl_int_sub(bmap->div[div][1],
3120                                         bmap->div[div][1], bmap->div[div][0]);
3121                         isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
3122                         neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
3123                         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
3124                         isl_int_add(bmap->div[div][1],
3125                                         bmap->div[div][1], bmap->div[div][0]);
3126                         if (!neg)
3127                                 return 0;
3128                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
3129                                                     bmap->n_div-div-1) != -1)
3130                                 return 0;
3131                 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
3132                         if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
3133                                 return 0;
3134                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
3135                                                     bmap->n_div-div-1) != -1)
3136                                 return 0;
3137                 } else
3138                         return 0;
3139         }
3140
3141         for (i = 0; i < bmap->n_div; ++i)
3142                 if (!isl_int_is_zero(bmap->div[i][1+pos]))
3143                         return 0;
3144
3145         return 1;
3146 }
3147
3148 /*
3149  * Remove divs that don't occur in any of the constraints or other divs.
3150  * These can arise when dropping some of the variables in a quast
3151  * returned by piplib.
3152  */
3153 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
3154 {
3155         int i;
3156
3157         if (!bmap)
3158                 return NULL;
3159
3160         for (i = bmap->n_div-1; i >= 0; --i) {
3161                 if (!div_is_redundant(bmap, i))
3162                         continue;
3163                 bmap = isl_basic_map_drop_div(bmap, i);
3164         }
3165         return bmap;
3166 }
3167
3168 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
3169 {
3170         bmap = remove_redundant_divs(bmap);
3171         if (!bmap)
3172                 return NULL;
3173         F_SET(bmap, ISL_BASIC_SET_FINAL);
3174         return bmap;
3175 }
3176
3177 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
3178 {
3179         return (struct isl_basic_set *)
3180                 isl_basic_map_finalize((struct isl_basic_map *)bset);
3181 }
3182
3183 struct isl_set *isl_set_finalize(struct isl_set *set)
3184 {
3185         int i;
3186
3187         if (!set)
3188                 return NULL;
3189         for (i = 0; i < set->n; ++i) {
3190                 set->p[i] = isl_basic_set_finalize(set->p[i]);
3191                 if (!set->p[i])
3192                         goto error;
3193         }
3194         return set;
3195 error:
3196         isl_set_free(set);
3197         return NULL;
3198 }
3199
3200 struct isl_map *isl_map_finalize(struct isl_map *map)
3201 {
3202         int i;
3203
3204         if (!map)
3205                 return NULL;
3206         for (i = 0; i < map->n; ++i) {
3207                 map->p[i] = isl_basic_map_finalize(map->p[i]);
3208                 if (!map->p[i])
3209                         goto error;
3210         }
3211         F_CLR(map, ISL_MAP_NORMALIZED);
3212         return map;
3213 error:
3214         isl_map_free(map);
3215         return NULL;
3216 }
3217
3218 struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
3219                 unsigned nparam, unsigned dim)
3220 {
3221         struct isl_basic_map *bmap;
3222         int i;
3223
3224         bmap = isl_basic_map_alloc(ctx, nparam, dim, dim, 0, dim, 0);
3225         if (!bmap)
3226                 goto error;
3227
3228         for (i = 0; i < dim; ++i) {
3229                 int j = isl_basic_map_alloc_equality(bmap);
3230                 if (j < 0)
3231                         goto error;
3232                 isl_seq_clr(bmap->eq[j],
3233                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
3234                 isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
3235                 isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
3236         }
3237         return isl_basic_map_finalize(bmap);
3238 error:
3239         isl_basic_map_free(bmap);
3240         return NULL;
3241 }
3242
3243 struct isl_map *isl_map_identity(struct isl_ctx *ctx,
3244                 unsigned nparam, unsigned dim)
3245 {
3246         struct isl_map *map = isl_map_alloc(ctx, nparam, dim, dim, 1,
3247                                                 ISL_MAP_DISJOINT);
3248         if (!map)
3249                 goto error;
3250         map = isl_map_add(map,
3251                 isl_basic_map_identity(ctx, nparam, dim));
3252         return map;
3253 error:
3254         isl_map_free(map);
3255         return NULL;
3256 }
3257
3258 int isl_set_is_equal(struct isl_set *set1, struct isl_set *set2)
3259 {
3260         return isl_map_is_equal((struct isl_map *)set1, (struct isl_map *)set2);
3261 }
3262
3263 int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
3264 {
3265         return isl_map_is_subset(
3266                         (struct isl_map *)set1, (struct isl_map *)set2);
3267 }
3268
3269 int isl_basic_map_is_subset(
3270                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3271 {
3272         int is_subset;
3273         struct isl_map *map1;
3274         struct isl_map *map2;
3275
3276         if (!bmap1 || !bmap2)
3277                 return -1;
3278
3279         map1 = isl_map_from_basic_map(isl_basic_map_copy(bmap1));
3280         map2 = isl_map_from_basic_map(isl_basic_map_copy(bmap2));
3281
3282         is_subset = isl_map_is_subset(map1, map2);
3283
3284         isl_map_free(map1);
3285         isl_map_free(map2);
3286
3287         return is_subset;
3288 }
3289
3290 int isl_basic_map_is_equal(
3291                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3292 {
3293         int is_subset;
3294
3295         if (!bmap1 || !bmap2)
3296                 return -1;
3297         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3298         if (is_subset != 1)
3299                 return is_subset;
3300         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3301         return is_subset;
3302 }
3303
3304 int isl_basic_set_is_equal(
3305                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3306 {
3307         return isl_basic_map_is_equal(
3308                 (struct isl_basic_map *)bset1, (struct isl_basic_map *)bset2);
3309 }
3310
3311 int isl_map_is_empty(struct isl_map *map)
3312 {
3313         int i;
3314         int is_empty;
3315
3316         if (!map)
3317                 return -1;
3318         for (i = 0; i < map->n; ++i) {
3319                 is_empty = isl_basic_map_is_empty(map->p[i]);
3320                 if (is_empty < 0)
3321                         return -1;
3322                 if (!is_empty)
3323                         return 0;
3324         }
3325         return 1;
3326 }
3327
3328 int isl_set_is_empty(struct isl_set *set)
3329 {
3330         return isl_map_is_empty((struct isl_map *)set);
3331 }
3332
3333 int isl_map_is_subset(struct isl_map *map1, struct isl_map *map2)
3334 {
3335         int i;
3336         int is_subset = 0;
3337         struct isl_map *diff;
3338
3339         if (!map1 || !map2)
3340                 return -1;
3341
3342         if (isl_map_is_empty(map1))
3343                 return 1;
3344
3345         if (isl_map_is_empty(map2))
3346                 return 0;
3347
3348         diff = isl_map_subtract(isl_map_copy(map1), isl_map_copy(map2));
3349         if (!diff)
3350                 return -1;
3351
3352         is_subset = isl_map_is_empty(diff);
3353         isl_map_free(diff);
3354
3355         return is_subset;
3356 }
3357
3358 int isl_map_is_equal(struct isl_map *map1, struct isl_map *map2)
3359 {
3360         int is_subset;
3361
3362         if (!map1 || !map2)
3363                 return -1;
3364         is_subset = isl_map_is_subset(map1, map2);
3365         if (is_subset != 1)
3366                 return is_subset;
3367         is_subset = isl_map_is_subset(map2, map1);
3368         return is_subset;
3369 }
3370
3371 int isl_basic_map_is_strict_subset(
3372                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3373 {
3374         int is_subset;
3375
3376         if (!bmap1 || !bmap2)
3377                 return -1;
3378         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3379         if (is_subset != 1)
3380                 return is_subset;
3381         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3382         if (is_subset == -1)
3383                 return is_subset;
3384         return !is_subset;
3385 }
3386
3387 static int basic_map_contains(struct isl_basic_map *bmap, struct isl_vec *vec)
3388 {
3389         int i;
3390         unsigned total;
3391         isl_int s;
3392
3393         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3394         if (total != vec->size)
3395                 return -1;
3396
3397         isl_int_init(s);
3398
3399         for (i = 0; i < bmap->n_eq; ++i) {
3400                 isl_seq_inner_product(vec->block.data, bmap->eq[i], total, &s);
3401                 if (!isl_int_is_zero(s)) {
3402                         isl_int_clear(s);
3403                         return 0;
3404                 }
3405         }
3406
3407         for (i = 0; i < bmap->n_ineq; ++i) {
3408                 isl_seq_inner_product(vec->block.data, bmap->ineq[i], total, &s);
3409                 if (isl_int_is_neg(s)) {
3410                         isl_int_clear(s);
3411                         return 0;
3412                 }
3413         }
3414
3415         isl_int_clear(s);
3416
3417         return 1;
3418 }
3419
3420 int isl_basic_map_is_universe(struct isl_basic_map *bmap)
3421 {
3422         if (!bmap)
3423                 return -1;
3424         return bmap->n_eq == 0 && bmap->n_ineq == 0;
3425 }
3426
3427 int isl_basic_map_is_empty(struct isl_basic_map *bmap)
3428 {
3429         struct isl_basic_set *bset = NULL;
3430         struct isl_vec *sample = NULL;
3431         int empty;
3432         unsigned total;
3433
3434         if (!bmap)
3435                 return -1;
3436
3437         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
3438                 return 1;
3439
3440         if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
3441                 struct isl_basic_map *copy = isl_basic_map_copy(bmap);
3442                 copy = isl_basic_map_convex_hull(copy);
3443                 empty = F_ISSET(copy, ISL_BASIC_MAP_EMPTY);
3444                 isl_basic_map_free(copy);
3445                 return empty;
3446         }
3447
3448         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3449         if (bmap->sample && bmap->sample->size == total) {
3450                 int contains = basic_map_contains(bmap, bmap->sample);
3451                 if (contains < 0)
3452                         return -1;
3453                 if (contains)
3454                         return 0;
3455         }
3456         bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
3457         if (!bset)
3458                 return -1;
3459         sample = isl_basic_set_sample(bset);
3460         if (!sample)
3461                 return -1;
3462         empty = sample->size == 0;
3463         if (bmap->sample)
3464                 isl_vec_free(bmap->ctx, bmap->sample);
3465         bmap->sample = sample;
3466
3467         return empty;
3468 }
3469
3470 struct isl_map *isl_basic_map_union(
3471         struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3472 {
3473         struct isl_map *map;
3474         if (!bmap1 || !bmap2)
3475                 return NULL;
3476
3477         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
3478         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
3479         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
3480
3481         map = isl_map_alloc(bmap1->ctx, bmap1->nparam,
3482                                 bmap1->n_in, bmap1->n_out, 2, 0);
3483         if (!map)
3484                 goto error;
3485         map = isl_map_add(map, bmap1);
3486         map = isl_map_add(map, bmap2);
3487         return map;
3488 error:
3489         isl_basic_map_free(bmap1);
3490         isl_basic_map_free(bmap2);
3491         return NULL;
3492 }
3493
3494 struct isl_set *isl_basic_set_union(
3495                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3496 {
3497         return (struct isl_set *)isl_basic_map_union(
3498                                             (struct isl_basic_map *)bset1,
3499                                             (struct isl_basic_map *)bset2);
3500 }
3501
3502 /* Order divs such that any div only depends on previous divs */
3503 static struct isl_basic_map *order_divs(struct isl_basic_map *bmap)
3504 {
3505         int i;
3506         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
3507
3508         for (i = 0; i < bmap->n_div; ++i) {
3509                 int pos;
3510                 pos = isl_seq_first_non_zero(bmap->div[i]+1+1+off+i,
3511                                                             bmap->n_div-i);
3512                 if (pos == -1)
3513                         continue;
3514                 swap_div(bmap, i, pos);
3515                 --i;
3516         }
3517         return bmap;
3518 }
3519
3520 static int find_div(struct isl_basic_map *dst,
3521                         struct isl_basic_map *src, unsigned div)
3522 {
3523         int i;
3524
3525         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3526
3527         for (i = 0; i < dst->n_div; ++i)
3528                 if (isl_seq_eq(dst->div[i], src->div[div], 1+1+total) &&
3529                     isl_seq_first_non_zero(dst->div[i]+1+1+total,
3530                                                 dst->n_div - src->n_div) == -1)
3531                         return i;
3532         return -1;
3533 }
3534
3535 struct isl_basic_map *isl_basic_map_align_divs(
3536                 struct isl_basic_map *dst, struct isl_basic_map *src)
3537 {
3538         int i;
3539         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3540
3541         if (src->n_div == 0)
3542                 return dst;
3543
3544         src = order_divs(src);
3545         dst = isl_basic_map_extend(dst, dst->nparam, dst->n_in,
3546                         dst->n_out, src->n_div, 0, 2 * src->n_div);
3547         if (!dst)
3548                 return NULL;
3549         for (i = 0; i < src->n_div; ++i) {
3550                 int j = find_div(dst, src, i);
3551                 if (j < 0) {
3552                         j = isl_basic_map_alloc_div(dst);
3553                         if (j < 0)
3554                                 goto error;
3555                         isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
3556                         isl_seq_clr(dst->div[j]+1+1+total,
3557                                                     dst->extra - src->n_div);
3558                         if (add_div_constraints(dst, j) < 0)
3559                                 goto error;
3560                 }
3561                 if (j != i)
3562                         swap_div(dst, i, j);
3563         }
3564         return dst;
3565 error:
3566         isl_basic_map_free(dst);
3567         return NULL;
3568 }
3569
3570 struct isl_map *isl_map_align_divs(struct isl_map *map)
3571 {
3572         int i;
3573
3574         map = isl_map_compute_divs(map);
3575         map = isl_map_cow(map);
3576         if (!map)
3577                 return NULL;
3578
3579         for (i = 1; i < map->n; ++i)
3580                 map->p[0] = isl_basic_map_align_divs(map->p[0], map->p[i]);
3581         for (i = 1; i < map->n; ++i)
3582                 map->p[i] = isl_basic_map_align_divs(map->p[i], map->p[0]);
3583
3584         F_CLR(map, ISL_MAP_NORMALIZED);
3585         return map;
3586 }
3587
3588 static struct isl_map *add_cut_constraint(struct isl_map *dst,
3589                 struct isl_basic_map *src, isl_int *c,
3590                 unsigned len, int oppose)
3591 {
3592         struct isl_basic_map *copy = NULL;
3593         int is_empty;
3594         int k;
3595         unsigned total;
3596
3597         copy = isl_basic_map_copy(src);
3598         copy = isl_basic_map_cow(copy);
3599         if (!copy)
3600                 goto error;
3601         copy = isl_basic_map_extend(copy,
3602                 copy->nparam, copy->n_in, copy->n_out, 0, 0, 1);
3603         k = isl_basic_map_alloc_inequality(copy);
3604         if (k < 0)
3605                 goto error;
3606         if (oppose)
3607                 isl_seq_neg(copy->ineq[k], c, len);
3608         else
3609                 isl_seq_cpy(copy->ineq[k], c, len);
3610         total = 1 + copy->nparam + copy->n_in + copy->n_out + copy->extra;
3611         isl_seq_clr(copy->ineq[k]+len, total - len);
3612         isl_inequality_negate(copy, k);
3613         copy = isl_basic_map_simplify(copy);
3614         copy = isl_basic_map_finalize(copy);
3615         is_empty = isl_basic_map_is_empty(copy);
3616         if (is_empty < 0)
3617                 goto error;
3618         if (!is_empty)
3619                 dst = isl_map_add(dst, copy);
3620         else
3621                 isl_basic_map_free(copy);
3622         return dst;
3623 error:
3624         isl_basic_map_free(copy);
3625         isl_map_free(dst);
3626         return NULL;
3627 }
3628
3629 static struct isl_map *subtract(struct isl_map *map, struct isl_basic_map *bmap)
3630 {
3631         int i, j, k;
3632         unsigned flags = 0;
3633         struct isl_map *rest = NULL;
3634         unsigned max;
3635         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3636
3637         assert(bmap);
3638
3639         if (!map)
3640                 goto error;
3641
3642         if (F_ISSET(map, ISL_MAP_DISJOINT))
3643                 FL_SET(flags, ISL_MAP_DISJOINT);
3644
3645         max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
3646         rest = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
3647                                 max, flags);
3648         if (!rest)
3649                 goto error;
3650
3651         for (i = 0; i < map->n; ++i) {
3652                 map->p[i] = isl_basic_map_align_divs(map->p[i], bmap);
3653                 if (!map->p[i])
3654                         goto error;
3655         }
3656
3657         for (j = 0; j < map->n; ++j)
3658                 map->p[j] = isl_basic_map_cow(map->p[j]);
3659
3660         for (i = 0; i < bmap->n_eq; ++i) {
3661                 for (j = 0; j < map->n; ++j) {
3662                         rest = add_cut_constraint(rest,
3663                                 map->p[j], bmap->eq[i], 1+total, 0);
3664                         if (!rest)
3665                                 goto error;
3666
3667                         rest = add_cut_constraint(rest,
3668                                 map->p[j], bmap->eq[i], 1+total, 1);
3669                         if (!rest)
3670                                 goto error;
3671
3672                         map->p[j] = isl_basic_map_extend(map->p[j],
3673                                 map->p[j]->nparam, map->p[j]->n_in,
3674                                 map->p[j]->n_out, 0, 1, 0);
3675                         if (!map->p[j])
3676                                 goto error;
3677                         k = isl_basic_map_alloc_equality(map->p[j]);
3678                         if (k < 0)
3679                                 goto error;
3680                         isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
3681                         isl_seq_clr(map->p[j]->eq[k]+1+total,
3682                                         map->p[j]->extra - bmap->n_div);
3683                 }
3684         }
3685
3686         for (i = 0; i < bmap->n_ineq; ++i) {
3687                 for (j = 0; j < map->n; ++j) {
3688                         rest = add_cut_constraint(rest,
3689                                 map->p[j], bmap->ineq[i], 1+total, 0);
3690                         if (!rest)
3691                                 goto error;
3692
3693                         map->p[j] = isl_basic_map_extend(map->p[j],
3694                                 map->p[j]->nparam, map->p[j]->n_in,
3695                                 map->p[j]->n_out, 0, 0, 1);
3696                         if (!map->p[j])
3697                                 goto error;
3698                         k = isl_basic_map_alloc_inequality(map->p[j]);
3699                         if (k < 0)
3700                                 goto error;
3701                         isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
3702                         isl_seq_clr(map->p[j]->ineq[k]+1+total,
3703                                         map->p[j]->extra - bmap->n_div);
3704                 }
3705         }
3706
3707         isl_map_free(map);
3708         return rest;
3709 error:
3710         isl_map_free(map);
3711         isl_map_free(rest);
3712         return NULL;
3713 }
3714
3715 struct isl_map *isl_map_subtract(struct isl_map *map1, struct isl_map *map2)
3716 {
3717         int i;
3718         if (!map1 || !map2)
3719                 goto error;
3720
3721         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
3722         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
3723         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
3724
3725         if (isl_map_is_empty(map2)) {
3726                 isl_map_free(map2);
3727                 return map1;
3728         }
3729
3730         map1 = isl_map_compute_divs(map1);
3731         map2 = isl_map_compute_divs(map2);
3732         if (!map1 || !map2)
3733                 goto error;
3734
3735         for (i = 0; map1 && i < map2->n; ++i)
3736                 map1 = subtract(map1, map2->p[i]);
3737
3738         isl_map_free(map2);
3739         return map1;
3740 error:
3741         isl_map_free(map1);
3742         isl_map_free(map2);
3743         return NULL;
3744 }
3745
3746 struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
3747 {
3748         return (struct isl_set *)
3749                 isl_map_subtract(
3750                         (struct isl_map *)set1, (struct isl_map *)set2);
3751 }
3752
3753 struct isl_set *isl_set_apply(struct isl_set *set, struct isl_map *map)
3754 {
3755         if (!set || !map)
3756                 goto error;
3757         isl_assert(set->ctx, set->dim == map->n_in, goto error);
3758         map = isl_map_intersect_domain(map, set);
3759         set = isl_map_range(map);
3760         return set;
3761 error:
3762         isl_set_free(set);
3763         isl_map_free(map);
3764         return NULL;
3765 }
3766
3767 /* There is no need to cow as removing empty parts doesn't change
3768  * the meaning of the set.
3769  */
3770 struct isl_map *isl_map_remove_empty_parts(struct isl_map *map)
3771 {
3772         int i;
3773
3774         if (!map)
3775                 return NULL;
3776
3777         for (i = map->n-1; i >= 0; --i) {
3778                 if (!F_ISSET(map->p[i], ISL_BASIC_MAP_EMPTY))
3779                         continue;
3780                 isl_basic_map_free(map->p[i]);
3781                 if (i != map->n-1) {
3782                         F_CLR(map, ISL_MAP_NORMALIZED);
3783                         map->p[i] = map->p[map->n-1];
3784                 }
3785                 map->n--;
3786         }
3787
3788         return map;
3789 }
3790
3791 struct isl_set *isl_set_remove_empty_parts(struct isl_set *set)
3792 {
3793         return (struct isl_set *)
3794                 isl_map_remove_empty_parts((struct isl_map *)set);
3795 }
3796
3797 struct isl_basic_set *isl_set_copy_basic_set(struct isl_set *set)
3798 {
3799         struct isl_basic_set *bset;
3800         if (!set || set->n == 0)
3801                 return NULL;
3802         bset = set->p[set->n-1];
3803         isl_assert(set->ctx, F_ISSET(bset, ISL_BASIC_SET_FINAL), return NULL);
3804         return isl_basic_set_copy(bset);
3805 }
3806
3807 struct isl_set *isl_set_drop_basic_set(struct isl_set *set,
3808                                                 struct isl_basic_set *bset)
3809 {
3810         int i;
3811
3812         if (!set || !bset)
3813                 goto error;
3814         for (i = set->n-1; i >= 0; --i) {
3815                 if (set->p[i] != bset)
3816                         continue;
3817                 set = isl_set_cow(set);
3818                 if (!set)
3819                         goto error;
3820                 isl_basic_set_free(set->p[i]);
3821                 if (i != set->n-1) {
3822                         F_CLR(set, ISL_SET_NORMALIZED);
3823                         set->p[i] = set->p[set->n-1];
3824                 }
3825                 set->n--;
3826                 return set;
3827         }
3828         isl_basic_set_free(bset);
3829         return set;
3830 error:
3831         isl_set_free(set);
3832         isl_basic_set_free(bset);
3833         return NULL;
3834 }
3835
3836 /* Given two _disjoint_ basic sets bset1 and bset2, check whether
3837  * for any common value of the parameters and dimensions preceding dim
3838  * in both basic sets, the values of dimension pos in bset1 are
3839  * smaller or larger then those in bset2.
3840  *
3841  * Returns
3842  *       1 if bset1 follows bset2
3843  *      -1 if bset1 precedes bset2
3844  *       0 if bset1 and bset2 are incomparable
3845  *      -2 if some error occurred.
3846  */
3847 int isl_basic_set_compare_at(struct isl_basic_set *bset1,
3848         struct isl_basic_set *bset2, int pos)
3849 {
3850         struct isl_basic_map *bmap1 = NULL;
3851         struct isl_basic_map *bmap2 = NULL;
3852         struct isl_ctx *ctx;
3853         struct isl_vec *obj;
3854         unsigned total;
3855         isl_int num, den;
3856         enum isl_lp_result res;
3857         int cmp;
3858
3859         if (!bset1 || !bset2)
3860                 return -2;
3861
3862         bmap1 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset1),
3863                                                 pos, bset1->dim - pos);
3864         bmap2 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset2),
3865                                                 pos, bset2->dim - pos);
3866         if (!bmap1 || !bmap2)
3867                 goto error;
3868         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
3869                         bmap1->n_in, bmap1->n_out + bmap2->n_out,
3870                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
3871         if (!bmap1)
3872                 goto error;
3873         total = bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
3874         ctx = bmap1->ctx;
3875         obj = isl_vec_alloc(ctx, total);
3876         isl_seq_clr(obj->block.data, total);
3877         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in], 1);
3878         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in+
3879                                         bmap1->n_out-bmap2->n_out], -1);
3880         if (!obj)
3881                 goto error;
3882         bmap1 = add_constraints(bmap1, bmap2, 0, bmap1->n_out - bmap2->n_out);
3883         isl_int_init(num);
3884         isl_int_init(den);
3885         res = isl_solve_lp(bmap1, 0, obj->block.data, ctx->one, &num, &den);
3886         if (res == isl_lp_empty)
3887                 cmp = 0;
3888         else if (res == isl_lp_ok && isl_int_is_pos(num))
3889                 cmp = 1;
3890         else if ((res == isl_lp_ok && isl_int_is_neg(num)) ||
3891                   res == isl_lp_unbounded)
3892                 cmp = -1;
3893         else
3894                 cmp = -2;
3895         isl_int_clear(num);
3896         isl_int_clear(den);
3897         isl_basic_map_free(bmap1);
3898         isl_vec_free(ctx, obj);
3899         return cmp;
3900 error:
3901         isl_basic_map_free(bmap1);
3902         isl_basic_map_free(bmap2);
3903         return -2;
3904 }
3905
3906 static int isl_basic_map_fast_has_fixed_var(struct isl_basic_map *bmap,
3907         unsigned pos, isl_int *val)
3908 {
3909         int i;
3910         int d;
3911         unsigned total;
3912
3913         if (!bmap)
3914                 return -1;
3915         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3916         for (i = 0, d = total-1; i < bmap->n_eq && d+1 > pos; ++i) {
3917                 for (; d+1 > pos; --d)
3918                         if (!isl_int_is_zero(bmap->eq[i][1+d]))
3919                                 break;
3920                 if (d != pos)
3921                         continue;
3922                 if (isl_seq_first_non_zero(bmap->eq[i]+1, d) != -1)
3923                         return 0;
3924                 if (isl_seq_first_non_zero(bmap->eq[i]+1+d+1, total-d-1) != -1)
3925                         return 0;
3926                 if (!isl_int_is_one(bmap->eq[i][1+d]))
3927                         return 0;
3928                 if (val)
3929                         isl_int_neg(*val, bmap->eq[i][0]);
3930                 return 1;
3931         }
3932         return 0;
3933 }
3934
3935 static int isl_map_fast_has_fixed_var(struct isl_map *map,
3936         unsigned pos, isl_int *val)
3937 {
3938         int i;
3939         isl_int v;
3940         isl_int tmp;
3941         int fixed;
3942
3943         if (!map)
3944                 return -1;
3945         if (map->n == 0)
3946                 return 0;
3947         if (map->n == 1)
3948                 return isl_basic_map_fast_has_fixed_var(map->p[0], pos, val);
3949         isl_int_init(v);
3950         isl_int_init(tmp);
3951         fixed = isl_basic_map_fast_has_fixed_var(map->p[0], pos, &v);
3952         for (i = 1; fixed == 1 && i < map->n; ++i) {
3953                 fixed = isl_basic_map_fast_has_fixed_var(map->p[i], pos, &tmp);
3954                 if (fixed == 1 && isl_int_ne(tmp, v))
3955                         fixed = 0;
3956         }
3957         if (val)
3958                 isl_int_set(*val, v);
3959         isl_int_clear(tmp);
3960         isl_int_clear(v);
3961         return fixed;
3962 }
3963
3964 static int isl_set_fast_has_fixed_var(struct isl_set *set, unsigned pos,
3965         isl_int *val)
3966 {
3967         return isl_map_fast_has_fixed_var((struct isl_map *)set, pos, val);
3968 }
3969
3970 /* Check if dimension dim has fixed value and if so and if val is not NULL,
3971  * then return this fixed value in *val.
3972  */
3973 int isl_set_fast_dim_is_fixed(struct isl_set *set, unsigned dim, isl_int *val)
3974 {
3975         return isl_set_fast_has_fixed_var(set, set->nparam + dim, val);
3976 }
3977
3978 /* Check if input variable in has fixed value and if so and if val is not NULL,
3979  * then return this fixed value in *val.
3980  */
3981 int isl_map_fast_input_is_fixed(struct isl_map *map, unsigned in, isl_int *val)
3982 {
3983         return isl_map_fast_has_fixed_var(map, map->nparam + in, val);
3984 }
3985
3986 /* Check if dimension dim has an (obvious) fixed lower bound and if so
3987  * and if val is not NULL, then return this lower bound in *val.
3988  */
3989 int isl_basic_set_fast_dim_has_fixed_lower_bound(struct isl_basic_set *bset,
3990         unsigned dim, isl_int *val)
3991 {
3992         int i, i_eq = -1, i_ineq = -1;
3993         isl_int *c;
3994         unsigned total;
3995
3996         if (!bset)
3997                 return -1;
3998         total = bset->nparam + bset->dim + bset->n_div;
3999         for (i = 0; i < bset->n_eq; ++i) {
4000                 if (isl_int_is_zero(bset->eq[i][1+bset->nparam+dim]))
4001                         continue;
4002                 if (i_eq != -1)
4003                         return 0;
4004                 i_eq = i;
4005         }
4006         for (i = 0; i < bset->n_ineq; ++i) {
4007                 if (!isl_int_is_pos(bset->ineq[i][1+bset->nparam+dim]))
4008                         continue;
4009                 if (i_eq != -1 || i_ineq != -1)
4010                         return 0;
4011                 i_ineq = i;
4012         }
4013         if (i_eq == -1 && i_ineq == -1)
4014                 return 0;
4015         c = i_eq != -1 ? bset->eq[i_eq] : bset->ineq[i_ineq];
4016         /* The coefficient should always be one due to normalization. */
4017         if (!isl_int_is_one(c[1+bset->nparam+dim]))
4018                 return 0;
4019         if (isl_seq_first_non_zero(c+1, bset->nparam+dim) != -1)
4020                 return 0;
4021         if (isl_seq_first_non_zero(c+1+bset->nparam+dim+1,
4022                                         total - bset->nparam - dim - 1) != -1)
4023                 return 0;
4024         if (val)
4025                 isl_int_neg(*val, c[0]);
4026         return 1;
4027 }
4028
4029 int isl_set_fast_dim_has_fixed_lower_bound(struct isl_set *set,
4030         unsigned dim, isl_int *val)
4031 {
4032         int i;
4033         isl_int v;
4034         isl_int tmp;
4035         int fixed;
4036
4037         if (!set)
4038                 return -1;
4039         if (set->n == 0)
4040                 return 0;
4041         if (set->n == 1)
4042                 return isl_basic_set_fast_dim_has_fixed_lower_bound(set->p[0],
4043                                                                 dim, val);
4044         isl_int_init(v);
4045         isl_int_init(tmp);
4046         fixed = isl_basic_set_fast_dim_has_fixed_lower_bound(set->p[0],
4047                                                                 dim, &v);
4048         for (i = 1; fixed == 1 && i < set->n; ++i) {
4049                 fixed = isl_basic_set_fast_dim_has_fixed_lower_bound(set->p[i],
4050                                                                 dim, &tmp);
4051                 if (fixed == 1 && isl_int_ne(tmp, v))
4052                         fixed = 0;
4053         }
4054         if (val)
4055                 isl_int_set(*val, v);
4056         isl_int_clear(tmp);
4057         isl_int_clear(v);
4058         return fixed;
4059 }
4060
4061 /* Remove all information from bset that is redundant in the context
4062  * of context.  In particular, equalities that are linear combinations
4063  * of those in context are remobed.  Then the inequalities that are
4064  * redundant in the context of the equalities and inequalities of
4065  * context are removed.
4066  */
4067 static struct isl_basic_set *uset_gist(struct isl_basic_set *bset,
4068         struct isl_basic_set *context)
4069 {
4070         int i;
4071         isl_int opt;
4072         struct isl_basic_set *combined;
4073         struct isl_ctx *ctx;
4074         enum isl_lp_result res = isl_lp_ok;
4075
4076         if (!bset || !context)
4077                 goto error;
4078
4079         if (context->n_eq > 0) {
4080                 struct isl_mat *T;
4081                 struct isl_mat *T2;
4082                 struct isl_ctx *ctx = context->ctx;
4083                 context = isl_basic_set_remove_equalities(context, &T, &T2);
4084                 if (!context)
4085                         goto error;
4086                 bset = isl_basic_set_preimage(ctx, bset, T);
4087                 bset = uset_gist(bset, context);
4088                 bset = isl_basic_set_preimage(ctx, bset, T2);
4089                 return bset;
4090         }
4091         combined = isl_basic_set_extend(isl_basic_set_copy(bset),
4092                         0, bset->dim, 0, context->n_eq, context->n_ineq);
4093         context = set_add_constraints(combined, context, 0);
4094         if (!context)
4095                 goto error;
4096         ctx = context->ctx;
4097         isl_int_init(opt);
4098         for (i = bset->n_ineq-1; i >= 0; --i) {
4099                 set_swap_inequality(context, i, context->n_ineq-1);
4100                 context->n_ineq--;
4101                 res = isl_solve_lp((struct isl_basic_map *)context, 0,
4102                         context->ineq[context->n_ineq]+1, ctx->one, &opt, NULL);
4103                 context->n_ineq++;
4104                 set_swap_inequality(context, i, context->n_ineq-1);
4105                 if (res == isl_lp_unbounded)
4106                         continue;
4107                 if (res == isl_lp_error)
4108                         break;
4109                 if (res == isl_lp_empty) {
4110                         bset = isl_basic_set_set_to_empty(bset);
4111                         break;
4112                 }
4113                 isl_int_add(opt, opt, context->ineq[i][0]);
4114                 if (!isl_int_is_neg(opt)) {
4115                         isl_basic_set_drop_inequality(context, i);
4116                         isl_basic_set_drop_inequality(bset, i);
4117                 }
4118         }
4119         isl_int_clear(opt);
4120         if (res == isl_lp_error)
4121                 goto error;
4122         isl_basic_set_free(context);
4123         return bset;
4124 error:
4125         isl_basic_set_free(context);
4126         isl_basic_set_free(bset);
4127         return NULL;
4128 }
4129
4130 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
4131         struct isl_basic_map *context)
4132 {
4133         struct isl_basic_set *bset;
4134
4135         if (!bmap || !context)
4136                 goto error;
4137
4138         context = isl_basic_map_align_divs(context, bmap);
4139         bmap = isl_basic_map_align_divs(bmap, context);
4140
4141         bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
4142                          isl_basic_map_underlying_set(context));
4143
4144         return isl_basic_map_overlying_set(bset, bmap);
4145 error:
4146         isl_basic_map_free(bmap);
4147         isl_basic_map_free(context);
4148         return NULL;
4149 }
4150
4151 struct isl_map *isl_map_gist(struct isl_map *map, struct isl_basic_map *context)
4152 {
4153         int i;
4154
4155         map = isl_map_cow(map);
4156         if (!map || !context)
4157                 return NULL;
4158         isl_assert(map->ctx, map->nparam == context->nparam, goto error);
4159         isl_assert(map->ctx, map->n_in == context->n_in, goto error);
4160         isl_assert(map->ctx, map->n_out == context->n_out, goto error);
4161         for (i = 0; i < map->n; ++i)
4162                 context = isl_basic_map_align_divs(context, map->p[i]);
4163         for (i = 0; i < map->n; ++i) {
4164                 map->p[i] = isl_basic_map_gist(map->p[i],
4165                                                 isl_basic_map_copy(context));
4166                 if (!map->p[i])
4167                         goto error;
4168         }
4169         isl_basic_map_free(context);
4170         F_CLR(map, ISL_MAP_NORMALIZED);
4171         return map;
4172 error:
4173         isl_map_free(map);
4174         isl_basic_map_free(context);
4175         return NULL;
4176 }
4177
4178 struct isl_set *isl_set_gist(struct isl_set *set, struct isl_basic_set *context)
4179 {
4180         return (struct isl_set *)isl_map_gist((struct isl_map *)set,
4181                                         (struct isl_basic_map *)context);
4182 }
4183
4184 struct constraint {
4185         unsigned        size;
4186         isl_int         *c;
4187 };
4188
4189 static int qsort_constraint_cmp(const void *p1, const void *p2)
4190 {
4191         const struct constraint *c1 = (const struct constraint *)p1;
4192         const struct constraint *c2 = (const struct constraint *)p2;
4193         unsigned size = isl_min(c1->size, c2->size);
4194         return isl_seq_cmp(c1->c, c2->c, size);
4195 }
4196
4197 static struct isl_basic_map *isl_basic_map_sort_constraints(
4198         struct isl_basic_map *bmap)
4199 {
4200         int i;
4201         struct constraint *c;
4202         unsigned total;
4203
4204         if (!bmap)
4205                 return NULL;
4206         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
4207         c = isl_alloc_array(bmap->ctx, struct constraint, bmap->n_ineq);
4208         if (!c)
4209                 goto error;
4210         for (i = 0; i < bmap->n_ineq; ++i) {
4211                 c[i].size = total;
4212                 c[i].c = bmap->ineq[i];
4213         }
4214         qsort(c, bmap->n_ineq, sizeof(struct constraint), qsort_constraint_cmp);
4215         for (i = 0; i < bmap->n_ineq; ++i)
4216                 bmap->ineq[i] = c[i].c;
4217         free(c);
4218         return bmap;
4219 error:
4220         isl_basic_map_free(bmap);
4221         return NULL;
4222 }
4223
4224 struct isl_basic_map *isl_basic_map_normalize(struct isl_basic_map *bmap)
4225 {
4226         if (!bmap)
4227                 return NULL;
4228         if (F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED))
4229                 return bmap;
4230         bmap = isl_basic_map_convex_hull(bmap);
4231         bmap = isl_basic_map_sort_constraints(bmap);
4232         F_SET(bmap, ISL_BASIC_MAP_NORMALIZED);
4233         return bmap;
4234 }
4235
4236 static int isl_basic_map_fast_cmp(const struct isl_basic_map *bmap1,
4237         const struct isl_basic_map *bmap2)
4238 {
4239         int i, cmp;
4240         unsigned total;
4241
4242         if (bmap1 == bmap2)
4243                 return 0;
4244         if (bmap1->nparam != bmap2->nparam)
4245                 return bmap1->nparam - bmap2->nparam;
4246         if (bmap1->n_in != bmap2->n_in)
4247                 return bmap1->n_in - bmap2->n_in;
4248         if (bmap1->n_out != bmap2->n_out)
4249                 return bmap1->n_out - bmap2->n_out;
4250         if (F_ISSET(bmap1, ISL_BASIC_MAP_EMPTY) &&
4251             F_ISSET(bmap2, ISL_BASIC_MAP_EMPTY))
4252                 return 0;
4253         if (F_ISSET(bmap1, ISL_BASIC_MAP_EMPTY))
4254                 return 1;
4255         if (F_ISSET(bmap2, ISL_BASIC_MAP_EMPTY))
4256                 return -1;
4257         if (bmap1->n_eq != bmap2->n_eq)
4258                 return bmap1->n_eq - bmap2->n_eq;
4259         if (bmap1->n_ineq != bmap2->n_ineq)
4260                 return bmap1->n_ineq - bmap2->n_ineq;
4261         if (bmap1->n_div != bmap2->n_div)
4262                 return bmap1->n_div - bmap2->n_div;
4263         total = bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
4264         for (i = 0; i < bmap1->n_eq; ++i) {
4265                 cmp = isl_seq_cmp(bmap1->eq[i], bmap2->eq[i], 1+total);
4266                 if (cmp)
4267                         return cmp;
4268         }
4269         for (i = 0; i < bmap1->n_ineq; ++i) {
4270                 cmp = isl_seq_cmp(bmap1->ineq[i], bmap2->ineq[i], 1+total);
4271                 if (cmp)
4272                         return cmp;
4273         }
4274         for (i = 0; i < bmap1->n_div; ++i) {
4275                 cmp = isl_seq_cmp(bmap1->div[i], bmap2->div[i], 1+1+total);
4276                 if (cmp)
4277                         return cmp;
4278         }
4279         return 0;
4280 }
4281
4282 static int isl_basic_map_fast_is_equal(struct isl_basic_map *bmap1,
4283         struct isl_basic_map *bmap2)
4284 {
4285         return isl_basic_map_fast_cmp(bmap1, bmap2) == 0;
4286 }
4287
4288 static int qsort_bmap_cmp(const void *p1, const void *p2)
4289 {
4290         const struct isl_basic_map *bmap1 = *(const struct isl_basic_map **)p1;
4291         const struct isl_basic_map *bmap2 = *(const struct isl_basic_map **)p2;
4292
4293         return isl_basic_map_fast_cmp(bmap1, bmap2);
4294 }
4295
4296 /* We normalize in place, but if anything goes wrong we need
4297  * to return NULL, so we need to make sure we don't change the
4298  * meaning of any possible other copies of map.
4299  */
4300 struct isl_map *isl_map_normalize(struct isl_map *map)
4301 {
4302         int i, j;
4303         struct isl_basic_map *bmap;
4304
4305         if (!map)
4306                 return NULL;
4307         if (F_ISSET(map, ISL_MAP_NORMALIZED))
4308                 return map;
4309         for (i = 0; i < map->n; ++i) {
4310                 bmap = isl_basic_map_normalize(isl_basic_map_copy(map->p[i]));
4311                 if (!bmap)
4312                         goto error;
4313                 isl_basic_map_free(map->p[i]);
4314                 map->p[i] = bmap;
4315         }
4316         qsort(map->p, map->n, sizeof(struct isl_basic_map *), qsort_bmap_cmp);
4317         F_SET(map, ISL_MAP_NORMALIZED);
4318         map = isl_map_remove_empty_parts(map);
4319         if (!map)
4320                 return NULL;
4321         for (i = map->n - 1; i >= 1; --i) {
4322                 if (!isl_basic_map_fast_is_equal(map->p[i-1], map->p[i]))
4323                         continue;
4324                 isl_basic_map_free(map->p[i-1]);
4325                 for (j = i; j < map->n; ++j)
4326                         map->p[j-1] = map->p[j];
4327                 map->n--;
4328         }
4329         return map;
4330 error:
4331         isl_map_free(map);
4332         return NULL;
4333
4334 }
4335
4336 struct isl_set *isl_set_normalize(struct isl_set *set)
4337 {
4338         return (struct isl_set *)isl_map_normalize((struct isl_map *)set);
4339 }
4340
4341 int isl_map_fast_is_equal(struct isl_map *map1, struct isl_map *map2)
4342 {
4343         int i;
4344         int equal;
4345
4346         if (!map1 || !map2)
4347                 return -1;
4348
4349         if (map1 == map2)
4350                 return 1;
4351         if (map1->nparam != map2->nparam)
4352                 return 0;
4353         if (map1->n_in != map2->n_in)
4354                 return 0;
4355         if (map1->n_out != map2->n_out)
4356                 return 0;
4357
4358         map1 = isl_map_copy(map1);
4359         map2 = isl_map_copy(map2);
4360         map1 = isl_map_normalize(map1);
4361         map2 = isl_map_normalize(map2);
4362         if (!map1 || !map2)
4363                 goto error;
4364         equal = map1->n == map2->n;
4365         for (i = 0; equal && i < map1->n; ++i) {
4366                 equal = isl_basic_map_fast_is_equal(map1->p[i], map2->p[i]);
4367                 if (equal < 0)
4368                         goto error;
4369         }
4370         isl_map_free(map1);
4371         isl_map_free(map2);
4372         return equal;
4373 error:
4374         isl_map_free(map1);
4375         isl_map_free(map2);
4376         return -1;
4377 }
4378
4379 int isl_set_fast_is_equal(struct isl_set *set1, struct isl_set *set2)
4380 {
4381         return isl_map_fast_is_equal((struct isl_map *)set1,
4382                                                 (struct isl_map *)set2);
4383 }