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