isl_map.c

   1 #include <string.h>
   2 #include <strings.h>
   3 #include "isl_ctx.h"
   4 #include "isl_blk.h"
   5 #include "isl_seq.h"
   6 #include "isl_set.h"
   7 #include "isl_map.h"
   8 #include "isl_map_private.h"
   9 #include "isl_map_piplib.h"
  10 #include "isl_sample.h"
  11 #include "isl_vec.h"
  12
  13 static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
  14                 struct isl_basic_map *bmap,
  15                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
  16                 unsigned n_eq, unsigned n_ineq)
  17 {
  18         int i;
  19         size_t row_size = 1 + nparam + n_in + n_out + extra;
  20
  21         bmap->block = isl_blk_alloc(ctx, (n_eq + n_ineq) * row_size);
  22         if (!bmap->block.data) {
  23                 free(bmap);
  24                 return NULL;
  25         }
  26
  27         bmap->eq = isl_alloc_array(ctx, isl_int *, n_eq + n_ineq);
  28         if (!bmap->eq) {
  29                 isl_blk_free(ctx, bmap->block);
  30                 free(bmap);
  31                 return NULL;
  32         }
  33
  34         if (extra == 0) {
  35                 bmap->block2.size = 0;
  36                 bmap->block2.data = NULL;
  37                 bmap->div = NULL;
  38         } else {
  39                 bmap->block2 = isl_blk_alloc(ctx, extra * (1 + row_size));
  40                 if (!bmap->block2.data) {
  41                         free(bmap->eq);
  42                         isl_blk_free(ctx, bmap->block);
  43                         free(bmap);
  44                         return NULL;
  45                 }
  46
  47                 bmap->div = isl_alloc_array(ctx, isl_int *, extra);
  48                 if (!bmap->div) {
  49                         isl_blk_free(ctx, bmap->block2);
  50                         free(bmap->eq);
  51                         isl_blk_free(ctx, bmap->block);
  52                         free(bmap);
  53                         return NULL;
  54                 }
  55         }
  56
  57         for (i = 0; i < n_eq + n_ineq; ++i)
  58                 bmap->eq[i] = bmap->block.data + i * row_size;
  59
  60         for (i = 0; i < extra; ++i)
  61                 bmap->div[i] = bmap->block2.data + i * (1 + row_size);
  62
  63         bmap->ref = 1;
  64         bmap->flags = 0;
  65         bmap->c_size = n_eq + n_ineq;
  66         bmap->ineq = bmap->eq + n_eq;
  67         bmap->nparam = nparam;
  68         bmap->n_in = n_in;
  69         bmap->n_out = n_out;
  70         bmap->extra = extra;
  71         bmap->n_eq = 0;
  72         bmap->n_ineq = 0;
  73         bmap->n_div = 0;
  74
  75         return bmap;
  76 }
  77
  78 struct isl_basic_set *isl_basic_set_alloc(struct isl_ctx *ctx,
  79                 unsigned nparam, unsigned dim, unsigned extra,
  80                 unsigned n_eq, unsigned n_ineq)
  81 {
  82         struct isl_basic_map *bmap;
  83         bmap = isl_basic_map_alloc(ctx, nparam, 0, dim, extra, n_eq, n_ineq);
  84         return (struct isl_basic_set *)bmap;
  85 }
  86
  87 struct isl_basic_map *isl_basic_map_alloc(struct isl_ctx *ctx,
  88                 unsigned nparam, unsigned in, unsigned out, unsigned extra,
  89                 unsigned n_eq, unsigned n_ineq)
  90 {
  91         struct isl_basic_map *bmap;
  92
  93         bmap = isl_alloc_type(ctx, struct isl_basic_map);
  94         if (!bmap)
  95                 return NULL;
  96
  97         return basic_map_init(ctx, bmap,
  98                                nparam, in, out, extra, n_eq, n_ineq);
  99 }
 100
 101 static void dup_constraints(struct isl_ctx *ctx,
 102                 struct isl_basic_map *dst, struct isl_basic_map *src)
 103 {
 104         int i;
 105         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
 106
 107         for (i = 0; i < src->n_eq; ++i) {
 108                 int j = isl_basic_map_alloc_equality(ctx, dst);
 109                 isl_seq_cpy(dst->eq[j], src->eq[i], 1+total);
 110         }
 111
 112         for (i = 0; i < src->n_ineq; ++i) {
 113                 int j = isl_basic_map_alloc_inequality(ctx, dst);
 114                 isl_seq_cpy(dst->ineq[j], src->ineq[i], 1+total);
 115         }
 116
 117         for (i = 0; i < src->n_div; ++i) {
 118                 int j = isl_basic_map_alloc_div(ctx, dst);
 119                 isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
 120         }
 121         F_SET(dst, ISL_BASIC_SET_FINAL);
 122 }
 123
 124 struct isl_basic_map *isl_basic_map_dup(struct isl_ctx *ctx,
 125                                         struct isl_basic_map *bmap)
 126 {
 127         struct isl_basic_map *dup;
 128
 129         dup = isl_basic_map_alloc(ctx, bmap->nparam,
 130                         bmap->n_in, bmap->n_out,
 131                         bmap->n_div, bmap->n_eq, bmap->n_ineq);
 132         if (!dup)
 133                 return NULL;
 134         dup_constraints(ctx, dup, bmap);
 135         return dup;
 136 }
 137
 138 struct isl_basic_set *isl_basic_set_dup(struct isl_ctx *ctx,
 139                                         struct isl_basic_set *bset)
 140 {
 141         struct isl_basic_map *dup;
 142
 143         dup = isl_basic_map_dup(ctx, (struct isl_basic_map *)bset);
 144         return (struct isl_basic_set *)dup;
 145 }
 146
 147 struct isl_basic_set *isl_basic_set_copy(struct isl_ctx *ctx,
 148                                         struct isl_basic_set *bset)
 149 {
 150         if (!bset)
 151                 return NULL;
 152
 153         if (F_ISSET(bset, ISL_BASIC_SET_FINAL)) {
 154                 bset->ref++;
 155                 return bset;
 156         }
 157         return isl_basic_set_dup(ctx, bset);
 158 }
 159
 160 struct isl_set *isl_set_copy(struct isl_ctx *ctx, struct isl_set *set)
 161 {
 162         if (!set)
 163                 return NULL;
 164
 165         set->ref++;
 166         return set;
 167 }
 168
 169 struct isl_basic_map *isl_basic_map_copy(struct isl_ctx *ctx,
 170                                         struct isl_basic_map *bmap)
 171 {
 172         if (!bmap)
 173                 return NULL;
 174
 175         if (F_ISSET(bmap, ISL_BASIC_SET_FINAL)) {
 176                 bmap->ref++;
 177                 return bmap;
 178         }
 179         return isl_basic_map_dup(ctx, bmap);
 180 }
 181
 182 struct isl_map *isl_map_copy(struct isl_ctx *ctx, struct isl_map *map)
 183 {
 184         if (!map)
 185                 return NULL;
 186
 187         map->ref++;
 188         return map;
 189 }
 190
 191 void isl_basic_map_free(struct isl_ctx *ctx, struct isl_basic_map *bmap)
 192 {
 193         if (!bmap)
 194                 return;
 195
 196         if (--bmap->ref > 0)
 197                 return;
 198
 199         free(bmap->div);
 200         isl_blk_free(ctx, bmap->block2);
 201         free(bmap->eq);
 202         isl_blk_free(ctx, bmap->block);
 203         free(bmap);
 204 }
 205
 206 void isl_basic_set_free(struct isl_ctx *ctx, struct isl_basic_set *bset)
 207 {
 208         isl_basic_map_free(ctx, (struct isl_basic_map *)bset);
 209 }
 210
 211 int isl_basic_map_alloc_equality(struct isl_ctx *ctx,
 212                 struct isl_basic_map *bmap)
 213 {
 214         if (!bmap)
 215                 return -1;
 216         isl_assert(ctx, bmap->n_eq + bmap->n_ineq < bmap->c_size, return -1);
 217         isl_assert(ctx, bmap->eq + bmap->n_eq <= bmap->ineq, return -1);
 218         if (bmap->eq + bmap->n_eq == bmap->ineq) {
 219                 isl_int *t;
 220                 int j = isl_basic_map_alloc_inequality(ctx, bmap);
 221                 if (j < 0)
 222                         return -1;
 223                 t = bmap->ineq[0];
 224                 bmap->ineq[0] = bmap->ineq[j];
 225                 bmap->ineq[j] = t;
 226                 bmap->n_ineq--;
 227                 bmap->ineq++;
 228         }
 229         return bmap->n_eq++;
 230 }
 231
 232 int isl_basic_map_free_equality(struct isl_ctx *ctx,
 233                 struct isl_basic_map *bmap, unsigned n)
 234 {
 235         isl_assert(ctx, n <= bmap->n_eq, return -1);
 236         bmap->n_eq -= n;
 237         return 0;
 238 }
 239
 240 int isl_basic_map_drop_equality(struct isl_ctx *ctx,
 241                 struct isl_basic_map *bmap, unsigned pos)
 242 {
 243         isl_int *t;
 244         isl_assert(ctx, pos < bmap->n_eq, return -1);
 245
 246         if (pos != bmap->n_eq - 1) {
 247                 t = bmap->eq[pos];
 248                 bmap->eq[pos] = bmap->eq[bmap->n_eq - 1];
 249                 bmap->eq[bmap->n_eq - 1] = t;
 250         }
 251         bmap->n_eq--;
 252         return 0;
 253 }
 254
 255 void isl_basic_map_inequality_to_equality(struct isl_ctx *ctx,
 256                 struct isl_basic_map *bmap, unsigned pos)
 257 {
 258         isl_int *t;
 259
 260         t = bmap->ineq[pos];
 261         bmap->ineq[pos] = bmap->ineq[0];
 262         bmap->ineq[0] = bmap->eq[bmap->n_eq];
 263         bmap->eq[bmap->n_eq] = t;
 264         bmap->n_eq++;
 265         bmap->n_ineq--;
 266         bmap->ineq++;
 267 }
 268
 269 int isl_basic_map_alloc_inequality(struct isl_ctx *ctx,
 270                 struct isl_basic_map *bmap)
 271 {
 272         isl_assert(ctx, (bmap->ineq - bmap->eq) + bmap->n_ineq < bmap->c_size,
 273                         return -1);
 274         return bmap->n_ineq++;
 275 }
 276
 277 int isl_basic_map_free_inequality(struct isl_ctx *ctx,
 278                 struct isl_basic_map *bmap, unsigned n)
 279 {
 280         isl_assert(ctx, n <= bmap->n_ineq, return -1);
 281         bmap->n_ineq -= n;
 282         return 0;
 283 }
 284
 285 int isl_basic_map_drop_inequality(struct isl_ctx *ctx,
 286                 struct isl_basic_map *bmap, unsigned pos)
 287 {
 288         isl_int *t;
 289         isl_assert(ctx, pos < bmap->n_ineq, return -1);
 290
 291         if (pos != bmap->n_ineq - 1) {
 292                 t = bmap->ineq[pos];
 293                 bmap->ineq[pos] = bmap->ineq[bmap->n_ineq - 1];
 294                 bmap->ineq[bmap->n_ineq - 1] = t;
 295         }
 296         bmap->n_ineq--;
 297         return 0;
 298 }
 299
 300 int isl_basic_map_alloc_div(struct isl_ctx *ctx,
 301                 struct isl_basic_map *bmap)
 302 {
 303         isl_assert(ctx, bmap->n_div < bmap->extra, return -1);
 304         return bmap->n_div++;
 305 }
 306
 307 int isl_basic_map_free_div(struct isl_ctx *ctx,
 308                 struct isl_basic_map *bmap, unsigned n)
 309 {
 310         isl_assert(ctx, n <= bmap->n_div, return -1);
 311         bmap->n_div -= n;
 312         return 0;
 313 }
 314
 315 /* Copy constraint from src to dst, putting the vars of src at offset
 316  * dim_off in dst and the divs of src at offset div_off in dst.
 317  * If both sets are actually map, then dim_off applies to the input
 318  * variables.
 319  */
 320 static void copy_constraint(struct isl_basic_map *dst_map, isl_int *dst,
 321                             struct isl_basic_map *src_map, isl_int *src,
 322                             unsigned in_off, unsigned out_off, unsigned div_off)
 323 {
 324         isl_int_set(dst[0], src[0]);
 325         isl_seq_cpy(dst+1, src+1, isl_min(dst_map->nparam, src_map->nparam));
 326         isl_seq_cpy(dst+1+dst_map->nparam+in_off,
 327                     src+1+src_map->nparam,
 328                     isl_min(dst_map->n_in-in_off, src_map->n_in));
 329         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+out_off,
 330                     src+1+src_map->nparam+src_map->n_in,
 331                     isl_min(dst_map->n_out-out_off, src_map->n_out));
 332         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+div_off,
 333                     src+1+src_map->nparam+src_map->n_in+src_map->n_out,
 334                     isl_min(dst_map->extra-div_off, src_map->extra));
 335 }
 336
 337 static void copy_div(struct isl_basic_map *dst_map, isl_int *dst,
 338                      struct isl_basic_map *src_map, isl_int *src,
 339                      unsigned in_off, unsigned out_off, unsigned div_off)
 340 {
 341         isl_int_set(dst[0], src[0]);
 342         copy_constraint(dst_map, dst+1, src_map, src+1, in_off, out_off, div_off);
 343 }
 344
 345 static struct isl_basic_map *add_constraints(struct isl_ctx *ctx,
 346                 struct isl_basic_map *bmap1,
 347                 struct isl_basic_map *bmap2, unsigned i_pos, unsigned o_pos)
 348 {
 349         int i;
 350         unsigned div_off = bmap1->n_div;
 351
 352         for (i = 0; i < bmap2->n_eq; ++i) {
 353                 int i1 = isl_basic_map_alloc_equality(ctx, bmap1);
 354                 if (i1 < 0)
 355                         goto error;
 356                 copy_constraint(bmap1, bmap1->eq[i1], bmap2, bmap2->eq[i],
 357                                 i_pos, o_pos, div_off);
 358         }
 359
 360         for (i = 0; i < bmap2->n_ineq; ++i) {
 361                 int i1 = isl_basic_map_alloc_inequality(ctx, bmap1);
 362                 if (i1 < 0)
 363                         goto error;
 364                 copy_constraint(bmap1, bmap1->ineq[i1], bmap2, bmap2->ineq[i],
 365                                 i_pos, o_pos, div_off);
 366         }
 367
 368         for (i = 0; i < bmap2->n_div; ++i) {
 369                 int i1 = isl_basic_map_alloc_div(ctx, bmap1);
 370                 if (i1 < 0)
 371                         goto error;
 372                 copy_div(bmap1, bmap1->div[i1], bmap2, bmap2->div[i],
 373                          i_pos, o_pos, div_off);
 374         }
 375
 376         isl_basic_map_free(ctx, bmap2);
 377
 378         return bmap1;
 379
 380 error:
 381         isl_basic_map_free(ctx, bmap1);
 382         isl_basic_map_free(ctx, bmap2);
 383         return NULL;
 384 }
 385
 386 struct isl_basic_map *isl_basic_map_extend(struct isl_ctx *ctx,
 387                 struct isl_basic_map *base,
 388                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
 389                 unsigned n_eq, unsigned n_ineq)
 390 {
 391         struct isl_basic_map *ext;
 392         int dims_ok;
 393
 394         base = isl_basic_map_cow(ctx, base);
 395         if (!base)
 396                 goto error;
 397
 398         dims_ok = base && base->nparam == nparam &&
 399                   base->n_in == n_in && base->n_out == n_out &&
 400                   base->extra >= base->n_div + extra;
 401
 402         if (dims_ok && n_eq == 0 && n_ineq == 0)
 403                 return base;
 404
 405         if (base) {
 406                 isl_assert(ctx, base->nparam <= nparam, goto error);
 407                 isl_assert(ctx, base->n_in <= n_in, goto error);
 408                 isl_assert(ctx, base->n_out <= n_out, goto error);
 409                 extra += base->extra;
 410                 n_eq += base->n_eq;
 411                 n_ineq += base->n_ineq;
 412         }
 413
 414         ext = isl_basic_map_alloc(ctx, nparam, n_in, n_out,
 415                                         extra, n_eq, n_ineq);
 416         if (!base)
 417                 return ext;
 418         if (!ext)
 419                 goto error;
 420
 421         ext = add_constraints(ctx, ext, base, 0, 0);
 422
 423         return ext;
 424
 425 error:
 426         isl_basic_map_free(ctx, base);
 427         return NULL;
 428 }
 429
 430 struct isl_basic_set *isl_basic_set_extend(struct isl_ctx *ctx,
 431                 struct isl_basic_set *base,
 432                 unsigned nparam, unsigned dim, unsigned extra,
 433                 unsigned n_eq, unsigned n_ineq)
 434 {
 435         return (struct isl_basic_set *)
 436                 isl_basic_map_extend(ctx, (struct isl_basic_map *)base,
 437                                         nparam, 0, dim, extra, n_eq, n_ineq);
 438 }
 439
 440 struct isl_basic_set *isl_basic_set_cow(struct isl_ctx *ctx,
 441                 struct isl_basic_set *bset)
 442 {
 443         return (struct isl_basic_set *)
 444                 isl_basic_map_cow(ctx, (struct isl_basic_map *)bset);
 445 }
 446
 447 struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
 448                 struct isl_basic_map *bmap)
 449 {
 450         if (!bmap)
 451                 return NULL;
 452
 453         if (bmap->ref > 1) {
 454                 bmap->ref--;
 455                 bmap = isl_basic_map_dup(ctx, bmap);
 456         }
 457         F_CLR(bmap, ISL_BASIC_SET_FINAL);
 458         return bmap;
 459 }
 460
 461 static struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
 462 {
 463         if (set->ref == 1)
 464                 return set;
 465         set->ref--;
 466         return isl_set_dup(ctx, set);
 467 }
 468
 469 struct isl_map *isl_map_cow(struct isl_ctx *ctx, struct isl_map *map)
 470 {
 471         if (map->ref == 1)
 472                 return map;
 473         map->ref--;
 474         return isl_map_dup(ctx, map);
 475 }
 476
 477 static void swap_vars(struct isl_blk blk, isl_int *a,
 478                         unsigned a_len, unsigned b_len)
 479 {
 480         isl_seq_cpy(blk.data, a+a_len, b_len);
 481         isl_seq_cpy(blk.data+b_len, a, a_len);
 482         isl_seq_cpy(a, blk.data, b_len+a_len);
 483 }
 484
 485 struct isl_basic_set *isl_basic_set_swap_vars(struct isl_ctx *ctx,
 486                 struct isl_basic_set *bset, unsigned n)
 487 {
 488         int i;
 489         struct isl_blk blk;
 490
 491         if (!bset)
 492                 goto error;
 493
 494         isl_assert(ctx, n <= bset->dim, goto error);
 495
 496         if (n == bset->dim)
 497                 return bset;
 498
 499         bset = isl_basic_set_cow(ctx, bset);
 500         if (!bset)
 501                 return NULL;
 502
 503         blk = isl_blk_alloc(ctx, bset->dim);
 504         if (!blk.data)
 505                 goto error;
 506
 507         for (i = 0; i < bset->n_eq; ++i)
 508                 swap_vars(blk,
 509                           bset->eq[i]+1+bset->nparam, n, bset->dim - n);
 510
 511         for (i = 0; i < bset->n_ineq; ++i)
 512                 swap_vars(blk,
 513                           bset->ineq[i]+1+bset->nparam, n, bset->dim - n);
 514
 515         for (i = 0; i < bset->n_div; ++i)
 516                 swap_vars(blk,
 517                           bset->div[i]+1+1+bset->nparam, n, bset->dim - n);
 518
 519         isl_blk_free(ctx, blk);
 520
 521         return bset;
 522
 523 error:
 524         isl_basic_set_free(ctx, bset);
 525         return NULL;
 526 }
 527
 528 struct isl_set *isl_set_swap_vars(struct isl_ctx *ctx,
 529                 struct isl_set *set, unsigned n)
 530 {
 531         int i;
 532         set = isl_set_cow(ctx, set);
 533         if (!set)
 534                 return NULL;
 535
 536         for (i = 0; i < set->n; ++i) {
 537                 set->p[i] = isl_basic_set_swap_vars(ctx, set->p[i], n);
 538                 if (!set->p[i]) {
 539                         isl_set_free(ctx, set);
 540                         return NULL;
 541                 }
 542         }
 543         return set;
 544 }
 545
 546 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
 547 {
 548         isl_seq_cpy(c, c + n, rem);
 549         isl_seq_clr(c + rem, n);
 550 }
 551
 552 struct isl_basic_set *isl_basic_set_drop_vars(struct isl_ctx *ctx,
 553                 struct isl_basic_set *bset, unsigned first, unsigned n)
 554 {
 555         int i;
 556
 557         if (!bset)
 558                 goto error;
 559
 560         isl_assert(ctx, first + n <= bset->dim, goto error);
 561
 562         if (n == 0)
 563                 return bset;
 564
 565         bset = isl_basic_set_cow(ctx, bset);
 566         if (!bset)
 567                 return NULL;
 568
 569         for (i = 0; i < bset->n_eq; ++i)
 570                 constraint_drop_vars(bset->eq[i]+1+bset->nparam+first, n,
 571                                      (bset->dim-first-n)+bset->extra);
 572
 573         for (i = 0; i < bset->n_ineq; ++i)
 574                 constraint_drop_vars(bset->ineq[i]+1+bset->nparam+first, n,
 575                                      (bset->dim-first-n)+bset->extra);
 576
 577         for (i = 0; i < bset->n_div; ++i)
 578                 constraint_drop_vars(bset->div[i]+1+1+bset->nparam+first, n,
 579                                      (bset->dim-first-n)+bset->extra);
 580
 581         bset->dim -= n;
 582         bset->extra += n;
 583
 584         return isl_basic_set_finalize(ctx, bset);
 585 error:
 586         isl_basic_set_free(ctx, bset);
 587         return NULL;
 588 }
 589
 590 /*
 591  * We don't cow, as the div is assumed to be redundant.
 592  */
 593 static struct isl_basic_map *isl_basic_map_drop_div(struct isl_ctx *ctx,
 594                 struct isl_basic_map *bmap, unsigned div)
 595 {
 596         int i;
 597         unsigned pos;
 598
 599         if (!bmap)
 600                 goto error;
 601
 602         pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
 603
 604         isl_assert(ctx, div < bmap->n_div, goto error);
 605
 606         for (i = 0; i < bmap->n_eq; ++i)
 607                 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
 608
 609         for (i = 0; i < bmap->n_ineq; ++i) {
 610                 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
 611                         isl_basic_map_drop_inequality(ctx, bmap, i);
 612                         --i;
 613                         continue;
 614                 }
 615                 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
 616         }
 617
 618         for (i = 0; i < bmap->n_div; ++i)
 619                 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
 620
 621         if (div != bmap->n_div - 1) {
 622                 int j;
 623                 isl_int *t = bmap->div[div];
 624
 625                 for (j = div; j < bmap->n_div - 1; ++j)
 626                         bmap->div[j] = bmap->div[j+1];
 627
 628                 bmap->div[bmap->n_div - 1] = t;
 629         }
 630         isl_basic_map_free_div(ctx, bmap, 1);
 631
 632         return bmap;
 633 error:
 634         isl_basic_map_free(ctx, bmap);
 635         return NULL;
 636 }
 637
 638 static void eliminate_div(struct isl_ctx *ctx, struct isl_basic_map *bmap,
 639                                 isl_int *eq, unsigned div)
 640 {
 641         int i;
 642         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
 643         unsigned len;
 644         len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 645
 646         for (i = 0; i < bmap->n_eq; ++i)
 647                 if (bmap->eq[i] != eq)
 648                         isl_seq_elim(bmap->eq[i], eq, pos, len, NULL);
 649
 650         for (i = 0; i < bmap->n_ineq; ++i)
 651                 isl_seq_elim(bmap->ineq[i], eq, pos, len, NULL);
 652
 653         /* We need to be careful about circular definitions,
 654          * so for now we just remove the definitions of other divs that
 655          * depend on this div and (possibly) recompute them later.
 656          */
 657         for (i = 0; i < bmap->n_div; ++i)
 658                 if (!isl_int_is_zero(bmap->div[i][0]) &&
 659                     !isl_int_is_zero(bmap->div[i][1 + pos]))
 660                         isl_seq_clr(bmap->div[i], 1 + len);
 661
 662         isl_basic_map_drop_div(ctx, bmap, div);
 663 }
 664
 665 struct isl_basic_map *isl_basic_map_set_to_empty(
 666                 struct isl_ctx *ctx, struct isl_basic_map *bmap)
 667 {
 668         int i = 0;
 669         unsigned total;
 670         if (!bmap)
 671                 goto error;
 672         total = bmap->nparam + bmap->n_in + bmap->n_out;
 673         isl_basic_map_free_div(ctx, bmap, bmap->n_div);
 674         isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
 675         if (bmap->n_eq > 0)
 676                 isl_basic_map_free_equality(ctx, bmap, bmap->n_eq-1);
 677         else {
 678                 isl_basic_map_alloc_equality(ctx, bmap);
 679                 if (i < 0)
 680                         goto error;
 681         }
 682         isl_int_set_si(bmap->eq[i][0], 1);
 683         isl_seq_clr(bmap->eq[i]+1, total);
 684         F_SET(bmap, ISL_BASIC_MAP_EMPTY);
 685         return isl_basic_map_finalize(ctx, bmap);
 686 error:
 687         isl_basic_map_free(ctx, bmap);
 688         return NULL;
 689 }
 690
 691 struct isl_basic_set *isl_basic_set_set_to_empty(
 692                 struct isl_ctx *ctx, struct isl_basic_set *bset)
 693 {
 694         return (struct isl_basic_set *)
 695                 isl_basic_map_set_to_empty(ctx, (struct isl_basic_map *)bset);
 696 }
 697
 698 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
 699 {
 700         isl_int *t = bmap->eq[a];
 701         bmap->eq[a] = bmap->eq[b];
 702         bmap->eq[b] = t;
 703 }
 704
 705 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
 706 {
 707         isl_int *t = bmap->ineq[a];
 708         bmap->ineq[a] = bmap->ineq[b];
 709         bmap->ineq[b] = t;
 710 }
 711
 712 static void swap_div(struct isl_basic_map *bmap, int a, int b)
 713 {
 714         int i;
 715         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
 716         isl_int *t = bmap->div[a];
 717         bmap->div[a] = bmap->div[b];
 718         bmap->div[b] = t;
 719
 720         for (i = 0; i < bmap->n_eq; ++i)
 721                 isl_int_swap(bmap->eq[i][1+off+a], bmap->eq[i][1+off+b]);
 722
 723         for (i = 0; i < bmap->n_ineq; ++i)
 724                 isl_int_swap(bmap->ineq[i][1+off+a], bmap->ineq[i][1+off+b]);
 725
 726         for (i = 0; i < bmap->n_div; ++i)
 727                 isl_int_swap(bmap->div[i][1+1+off+a], bmap->div[i][1+1+off+b]);
 728 }
 729
 730 struct isl_basic_map *isl_basic_map_gauss(struct isl_ctx *ctx,
 731         struct isl_basic_map *bmap, int *progress)
 732 {
 733         int k;
 734         int done;
 735         int last_var;
 736         unsigned total_var;
 737         unsigned total;
 738
 739         if (!bmap)
 740                 return NULL;
 741
 742         total_var = bmap->nparam + bmap->n_in + bmap->n_out;
 743         total = total_var + bmap->n_div;
 744
 745         last_var = total - 1;
 746         for (done = 0; done < bmap->n_eq; ++done) {
 747                 for (; last_var >= 0; --last_var) {
 748                         for (k = done; k < bmap->n_eq; ++k)
 749                                 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
 750                                         break;
 751                         if (k < bmap->n_eq)
 752                                 break;
 753                 }
 754                 if (last_var < 0)
 755                         break;
 756                 if (k != done)
 757                         swap_equality(bmap, k, done);
 758                 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
 759                         isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
 760
 761                 for (k = 0; k < bmap->n_eq; ++k) {
 762                         if (k == done)
 763                                 continue;
 764                         if (isl_int_is_zero(bmap->eq[k][1+last_var]))
 765                                 continue;
 766                         if (progress)
 767                                 *progress = 1;
 768                         isl_seq_elim(bmap->eq[k], bmap->eq[done],
 769                                         1+last_var, 1+total, NULL);
 770                 }
 771
 772                 for (k = 0; k < bmap->n_ineq; ++k) {
 773                         if (isl_int_is_zero(bmap->ineq[k][1+last_var]))
 774                                 continue;
 775                         if (progress)
 776                                 *progress = 1;
 777                         isl_seq_elim(bmap->ineq[k], bmap->eq[done],
 778                                         1+last_var, 1+total, NULL);
 779                 }
 780
 781                 for (k = 0; k < bmap->n_div; ++k) {
 782                         if (isl_int_is_zero(bmap->div[k][0]))
 783                                 continue;
 784                         if (isl_int_is_zero(bmap->div[k][1+1+last_var]))
 785                                 continue;
 786                         if (progress)
 787                                 *progress = 1;
 788                         isl_seq_elim(bmap->div[k]+1, bmap->eq[done],
 789                                         1+last_var, 1+total, &bmap->div[k][0]);
 790                 }
 791
 792                 if (last_var >= total_var &&
 793                     isl_int_is_zero(bmap->div[last_var - total_var][0])) {
 794                         unsigned div = last_var - total_var;
 795                         isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
 796                         isl_int_set_si(bmap->div[div][1+1+last_var], 0);
 797                         isl_int_set(bmap->div[div][0],
 798                                     bmap->eq[done][1+last_var]);
 799                 }
 800         }
 801         if (done == bmap->n_eq)
 802                 return bmap;
 803         for (k = done; k < bmap->n_eq; ++k) {
 804                 if (isl_int_is_zero(bmap->eq[k][0]))
 805                         continue;
 806                 return isl_basic_map_set_to_empty(ctx, bmap);
 807         }
 808         isl_basic_map_free_equality(ctx, bmap, bmap->n_eq-done);
 809         return bmap;
 810 }
 811
 812 struct isl_basic_set *isl_basic_set_gauss(struct isl_ctx *ctx,
 813         struct isl_basic_set *bset, int *progress)
 814 {
 815         return (struct isl_basic_set*)isl_basic_map_gauss(ctx,
 816                         (struct isl_basic_map *)bset, progress);
 817 }
 818
 819 static unsigned int round_up(unsigned int v)
 820 {
 821         int old_v = v;
 822
 823         while (v) {
 824                 old_v = v;
 825                 v ^= v & -v;
 826         }
 827         return old_v << 1;
 828 }
 829
 830 static int hash_index(int *index, unsigned int size, int bits,
 831                         struct isl_basic_map *bmap, int k)
 832 {
 833         int h;
 834         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 835         u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
 836         for (h = hash; index[h]; h = (h+1) % size)
 837                 if (k != index[h]-1 &&
 838                     isl_seq_eq(bmap->ineq[k]+1,
 839                                bmap->ineq[index[h]-1]+1, total))
 840                         break;
 841         return h;
 842 }
 843
 844 static struct isl_basic_map *remove_duplicate_constraints(
 845         struct isl_ctx *ctx,
 846         struct isl_basic_map *bmap, int *progress)
 847 {
 848         unsigned int size;
 849         int *index;
 850         int k, l, h;
 851         int bits;
 852         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
 853         isl_int sum;
 854
 855         if (bmap->n_ineq <= 1)
 856                 return bmap;
 857
 858         size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
 859         bits = ffs(size) - 1;
 860         index = isl_alloc_array(ctx, int, size);
 861         memset(index, 0, size * sizeof(int));
 862         if (!index)
 863                 return bmap;
 864
 865         index[isl_seq_hash(bmap->ineq[0]+1, total, bits)] = 1;
 866         for (k = 1; k < bmap->n_ineq; ++k) {
 867                 h = hash_index(index, size, bits, bmap, k);
 868                 if (!index[h]) {
 869                         index[h] = k+1;
 870                         continue;
 871                 }
 872                 *progress = 1;
 873                 l = index[h] - 1;
 874                 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
 875                         swap_inequality(bmap, k, l);
 876                 isl_basic_map_drop_inequality(ctx, bmap, k);
 877                 --k;
 878         }
 879         isl_int_init(sum);
 880         for (k = 0; k < bmap->n_ineq-1; ++k) {
 881                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
 882                 h = hash_index(index, size, bits, bmap, k);
 883                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
 884                 if (!index[h])
 885                         continue;
 886                 l = index[h] - 1;
 887                 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
 888                 if (isl_int_is_pos(sum))
 889                         continue;
 890                 if (isl_int_is_zero(sum)) {
 891                         /* We need to break out of the loop after these
 892                          * changes since the contents of the hash
 893                          * will no longer be valid.
 894                          * Plus, we probably we want to regauss first.
 895                          */
 896                         isl_basic_map_drop_inequality(ctx, bmap, l);
 897                         isl_basic_map_inequality_to_equality(ctx, bmap, k);
 898                 } else
 899                         bmap = isl_basic_map_set_to_empty(ctx, bmap);
 900                 break;
 901         }
 902         isl_int_clear(sum);
 903
 904         free(index);
 905         return bmap;
 906 }
 907
 908 static struct isl_basic_map *remove_duplicate_divs(struct isl_ctx *ctx,
 909         struct isl_basic_map *bmap, int *progress)
 910 {
 911         unsigned int size;
 912         int *index;
 913         int k, l, h;
 914         int bits;
 915         struct isl_blk eq;
 916         unsigned total_var = bmap->nparam + bmap->n_in + bmap->n_out;
 917         unsigned total = total_var + bmap->n_div;
 918
 919         if (bmap->n_div <= 1)
 920                 return bmap;
 921
 922         for (k = bmap->n_div - 1; k >= 0; --k)
 923                 if (!isl_int_is_zero(bmap->div[k][0]))
 924                         break;
 925         if (k <= 0)
 926                 return bmap;
 927
 928         size = round_up(4 * bmap->n_div / 3 - 1);
 929         bits = ffs(size) - 1;
 930         index = isl_alloc_array(ctx, int, size);
 931         memset(index, 0, size * sizeof(int));
 932         if (!index)
 933                 return bmap;
 934         eq = isl_blk_alloc(ctx, 1+total);
 935         if (!eq.data)
 936                 goto out;
 937
 938         isl_seq_clr(eq.data, 1+total);
 939         index[isl_seq_hash(bmap->div[k], 2+total, bits)] = k + 1;
 940         for (--k; k >= 0; --k) {
 941                 u_int32_t hash;
 942
 943                 if (isl_int_is_zero(bmap->div[k][0]))
 944                         continue;
 945
 946                 hash = isl_seq_hash(bmap->div[k], 2+total, bits);
 947                 for (h = hash; index[h]; h = (h+1) % size)
 948                         if (isl_seq_eq(bmap->div[k],
 949                                        bmap->div[index[h]-1], 2+total))
 950                                 break;
 951                 if (index[h]) {
 952                         *progress = 1;
 953                         l = index[h] - 1;
 954                         isl_int_set_si(eq.data[1+total_var+k], -1);
 955                         isl_int_set_si(eq.data[1+total_var+l], 1);
 956                         eliminate_div(ctx, bmap, eq.data, l);
 957                         isl_int_set_si(eq.data[1+total_var+k], 0);
 958                         isl_int_set_si(eq.data[1+total_var+l], 0);
 959                 }
 960                 index[h] = k+1;
 961         }
 962
 963         isl_blk_free(ctx, eq);
 964 out:
 965         free(index);
 966         return bmap;
 967 }
 968
 969 static struct isl_basic_map *eliminate_divs(
 970                 struct isl_ctx *ctx, struct isl_basic_map *bmap,
 971                 int *progress)
 972 {
 973         int d;
 974         int i;
 975         int modified = 0;
 976         unsigned off;
 977
 978         if (!bmap)
 979                 return NULL;
 980
 981         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
 982
 983         for (d = bmap->n_div - 1; d >= 0 ; --d) {
 984                 for (i = 0; i < bmap->n_eq; ++i) {
 985                         if (!isl_int_is_one(bmap->eq[i][off + d]) &&
 986                             !isl_int_is_negone(bmap->eq[i][off + d]))
 987                                 continue;
 988                         modified = 1;
 989                         *progress = 1;
 990                         eliminate_div(ctx, bmap, bmap->eq[i], d);
 991                         isl_basic_map_drop_equality(ctx, bmap, i);
 992                 }
 993         }
 994         if (modified)
 995                 return eliminate_divs(ctx, bmap, progress);
 996         return bmap;
 997 }
 998
 999 static struct isl_basic_map *normalize_constraints(
1000                 struct isl_ctx *ctx, struct isl_basic_map *bmap)
1001 {
1002         int i;
1003         isl_int gcd;
1004         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1005
1006         isl_int_init(gcd);
1007         for (i = bmap->n_eq - 1; i >= 0; --i) {
1008                 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
1009                 if (isl_int_is_zero(gcd)) {
1010                         if (!isl_int_is_zero(bmap->eq[i][0])) {
1011                                 bmap = isl_basic_map_set_to_empty(ctx, bmap);
1012                                 break;
1013                         }
1014                         isl_basic_map_drop_equality(ctx, bmap, i);
1015                         continue;
1016                 }
1017                 if (isl_int_is_one(gcd))
1018                         continue;
1019                 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
1020                         bmap = isl_basic_map_set_to_empty(ctx, bmap);
1021                         break;
1022                 }
1023                 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
1024         }
1025
1026         for (i = bmap->n_ineq - 1; i >= 0; --i) {
1027                 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
1028                 if (isl_int_is_zero(gcd)) {
1029                         if (isl_int_is_neg(bmap->ineq[i][0])) {
1030                                 bmap = isl_basic_map_set_to_empty(ctx, bmap);
1031                                 break;
1032                         }
1033                         isl_basic_map_drop_inequality(ctx, bmap, i);
1034                         continue;
1035                 }
1036                 if (isl_int_is_one(gcd))
1037                         continue;
1038                 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
1039                 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
1040         }
1041         isl_int_clear(gcd);
1042
1043         return bmap;
1044 }
1045
1046 struct isl_basic_map *isl_basic_map_simplify(
1047                 struct isl_ctx *ctx, struct isl_basic_map *bmap)
1048 {
1049         int progress = 1;
1050         if (!bmap)
1051                 return NULL;
1052         while (progress) {
1053                 progress = 0;
1054                 bmap = normalize_constraints(ctx, bmap);
1055                 bmap = eliminate_divs(ctx, bmap, &progress);
1056                 bmap = isl_basic_map_gauss(ctx, bmap, &progress);
1057                 bmap = remove_duplicate_divs(ctx, bmap, &progress);
1058                 bmap = remove_duplicate_constraints(ctx, bmap, &progress);
1059         }
1060         return bmap;
1061 }
1062
1063 struct isl_basic_set *isl_basic_set_simplify(
1064                 struct isl_ctx *ctx, struct isl_basic_set *bset)
1065 {
1066         return (struct isl_basic_set *)
1067                 isl_basic_map_simplify(ctx,
1068                         (struct isl_basic_map *)bset);
1069 }
1070
1071 static void dump_term(struct isl_basic_map *bmap,
1072                         isl_int c, int pos, FILE *out)
1073 {
1074         unsigned in = bmap->n_in;
1075         unsigned dim = bmap->n_in + bmap->n_out;
1076         if (!pos)
1077                 isl_int_print(out, c);
1078         else {
1079                 if (!isl_int_is_one(c))
1080                         isl_int_print(out, c);
1081                 if (pos < 1 + bmap->nparam)
1082                         fprintf(out, "p%d", pos - 1);
1083                 else if (pos < 1 + bmap->nparam + in)
1084                         fprintf(out, "i%d", pos - 1 - bmap->nparam);
1085                 else if (pos < 1 + bmap->nparam + dim)
1086                         fprintf(out, "o%d", pos - 1 - bmap->nparam - in);
1087                 else
1088                         fprintf(out, "e%d", pos - 1 - bmap->nparam - dim);
1089         }
1090 }
1091
1092 static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
1093                                 int sign, FILE *out)
1094 {
1095         int i;
1096         int first;
1097         unsigned len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1098         isl_int v;
1099
1100         isl_int_init(v);
1101         for (i = 0, first = 1; i < len; ++i) {
1102                 if (isl_int_sgn(c[i]) * sign <= 0)
1103                         continue;
1104                 if (!first)
1105                         fprintf(out, " + ");
1106                 first = 0;
1107                 isl_int_abs(v, c[i]);
1108                 dump_term(bmap, v, i, out);
1109         }
1110         isl_int_clear(v);
1111         if (first)
1112                 fprintf(out, "0");
1113 }
1114
1115 static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
1116                                 const char *op, FILE *out, int indent)
1117 {
1118         fprintf(out, "%*s", indent, "");
1119
1120         dump_constraint_sign(bmap, c, 1, out);
1121         fprintf(out, " %s ", op);
1122         dump_constraint_sign(bmap, c, -1, out);
1123
1124         fprintf(out, "\n");
1125 }
1126
1127 static void dump_constraints(struct isl_basic_map *bmap,
1128                                 isl_int **c, unsigned n,
1129                                 const char *op, FILE *out, int indent)
1130 {
1131         int i;
1132
1133         for (i = 0; i < n; ++i)
1134                 dump_constraint(bmap, c[i], op, out, indent);
1135 }
1136
1137 static void dump_affine(struct isl_basic_map *bmap, isl_int *exp, FILE *out)
1138 {
1139         int j;
1140         int first = 1;
1141         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1142
1143         for (j = 0; j < 1 + total; ++j) {
1144                 if (isl_int_is_zero(exp[j]))
1145                         continue;
1146                 if (!first && isl_int_is_pos(exp[j]))
1147                         fprintf(out, "+");
1148                 dump_term(bmap, exp[j], j, out);
1149                 first = 0;
1150         }
1151 }
1152
1153 static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
1154 {
1155         int i;
1156
1157         dump_constraints(bmap, bmap->eq, bmap->n_eq, "=", out, indent);
1158         dump_constraints(bmap, bmap->ineq, bmap->n_ineq, ">=", out, indent);
1159
1160         for (i = 0; i < bmap->n_div; ++i) {
1161                 fprintf(out, "%*s", indent, "");
1162                 fprintf(out, "e%d = [(", i);
1163                 dump_affine(bmap, bmap->div[i]+1, out);
1164                 fprintf(out, ")/");
1165                 isl_int_print(out, bmap->div[i][0]);
1166                 fprintf(out, "]\n");
1167         }
1168 }
1169
1170 void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
1171                                 FILE *out, int indent)
1172 {
1173         fprintf(out, "%*s", indent, "");
1174         fprintf(out, "nparam: %d, dim: %d, extra: %d\n",
1175                         bset->nparam, bset->dim, bset->extra);
1176         dump((struct isl_basic_map *)bset, out, indent);
1177 }
1178
1179 void isl_basic_map_dump(struct isl_ctx *ctx, struct isl_basic_map *bmap,
1180                                 FILE *out, int indent)
1181 {
1182         fprintf(out, "%*s", indent, "");
1183         fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d\n",
1184                 bmap->ref,
1185                 bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra);
1186         dump(bmap, out, indent);
1187 }
1188
1189 int isl_inequality_negate(struct isl_ctx *ctx,
1190                 struct isl_basic_map *bmap, unsigned pos)
1191 {
1192         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1193         isl_assert(ctx, pos < bmap->n_ineq, return -1);
1194         isl_seq_neg(bmap->ineq[pos], bmap->ineq[pos], 1 + total);
1195         isl_int_sub_ui(bmap->ineq[pos][0], bmap->ineq[pos][0], 1);
1196         return 0;
1197 }
1198
1199 struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
1200                 unsigned nparam, unsigned dim, int n, unsigned flags)
1201 {
1202         struct isl_set *set;
1203
1204         set = isl_alloc(ctx, struct isl_set,
1205                         sizeof(struct isl_set) +
1206                         n * sizeof(struct isl_basic_set *));
1207         if (!set)
1208                 return NULL;
1209
1210         set->ref = 1;
1211         set->size = n;
1212         set->n = 0;
1213         set->nparam = nparam;
1214         set->zero = 0;
1215         set->dim = dim;
1216         set->flags = flags;
1217         return set;
1218 }
1219
1220 struct isl_set *isl_set_dup(struct isl_ctx *ctx, struct isl_set *set)
1221 {
1222         int i;
1223         struct isl_set *dup;
1224
1225         dup = isl_set_alloc(ctx, set->nparam, set->dim, set->n, set->flags);
1226         if (!dup)
1227                 return NULL;
1228         for (i = 0; i < set->n; ++i)
1229                 dup = isl_set_add(ctx, dup,
1230                                   isl_basic_set_dup(ctx, set->p[i]));
1231         return dup;
1232 }
1233
1234 struct isl_set *isl_set_from_basic_set(struct isl_ctx *ctx,
1235                                 struct isl_basic_set *bset)
1236 {
1237         struct isl_set *set;
1238
1239         set = isl_set_alloc(ctx, bset->nparam, bset->dim, 1, ISL_MAP_DISJOINT);
1240         if (!set) {
1241                 isl_basic_set_free(ctx, bset);
1242                 return NULL;
1243         }
1244         return isl_set_add(ctx, set, bset);
1245 }
1246
1247 struct isl_map *isl_map_from_basic_map(struct isl_ctx *ctx,
1248                                 struct isl_basic_map *bmap)
1249 {
1250         struct isl_map *map;
1251
1252         map = isl_map_alloc(ctx, bmap->nparam, bmap->n_in, bmap->n_out, 1,
1253                                 ISL_MAP_DISJOINT);
1254         if (!map) {
1255                 isl_basic_map_free(ctx, bmap);
1256                 return NULL;
1257         }
1258         return isl_map_add(ctx, map, bmap);
1259 }
1260
1261 struct isl_set *isl_set_add(struct isl_ctx *ctx, struct isl_set *set,
1262                                 struct isl_basic_set *bset)
1263 {
1264         if (!bset || !set)
1265                 goto error;
1266         isl_assert(ctx, set->nparam == bset->nparam, goto error);
1267         isl_assert(ctx, set->dim == bset->dim, goto error);
1268         isl_assert(ctx, set->n < set->size, goto error);
1269         set->p[set->n] = bset;
1270         set->n++;
1271         return set;
1272 error:
1273         if (set)
1274                 isl_set_free(ctx, set);
1275         if (bset)
1276                 isl_basic_set_free(ctx, bset);
1277         return NULL;
1278 }
1279
1280 void isl_set_free(struct isl_ctx *ctx, struct isl_set *set)
1281 {
1282         int i;
1283
1284         if (!set)
1285                 return;
1286
1287         if (--set->ref > 0)
1288                 return;
1289
1290         for (i = 0; i < set->n; ++i)
1291                 isl_basic_set_free(ctx, set->p[i]);
1292         free(set);
1293 }
1294
1295 void isl_set_dump(struct isl_ctx *ctx, struct isl_set *set, FILE *out,
1296                   int indent)
1297 {
1298         int i;
1299
1300         if (!set) {
1301                 fprintf(out, "null set\n");
1302                 return;
1303         }
1304
1305         fprintf(out, "%*s", indent, "");
1306         fprintf(out, "ref: %d, n: %d\n", set->ref, set->n);
1307         for (i = 0; i < set->n; ++i) {
1308                 fprintf(out, "%*s", indent, "");
1309                 fprintf(out, "basic set %d:\n", i);
1310                 isl_basic_set_dump(ctx, set->p[i], out, indent+4);
1311         }
1312 }
1313
1314 void isl_map_dump(struct isl_ctx *ctx, struct isl_map *map, FILE *out,
1315                   int indent)
1316 {
1317         int i;
1318
1319         fprintf(out, "%*s", indent, "");
1320         fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d\n",
1321                         map->ref, map->n, map->nparam, map->n_in, map->n_out);
1322         for (i = 0; i < map->n; ++i) {
1323                 fprintf(out, "%*s", indent, "");
1324                 fprintf(out, "basic map %d:\n", i);
1325                 isl_basic_map_dump(ctx, map->p[i], out, indent+4);
1326         }
1327 }
1328
1329 struct isl_basic_map *isl_basic_map_intersect_domain(
1330                 struct isl_ctx *ctx, struct isl_basic_map *bmap,
1331                 struct isl_basic_set *bset)
1332 {
1333         struct isl_basic_map *bmap_domain;
1334
1335         if (!bmap || !bset)
1336                 goto error;
1337
1338         isl_assert(ctx, bset->dim == bmap->n_in, goto error);
1339         isl_assert(ctx, bset->nparam == bmap->nparam, goto error);
1340
1341         bmap = isl_basic_map_extend(ctx, bmap,
1342                         bset->nparam, bmap->n_in, bmap->n_out,
1343                         bset->n_div, bset->n_eq, bset->n_ineq);
1344         if (!bmap)
1345                 goto error;
1346         bmap_domain = isl_basic_map_from_basic_set(ctx, bset,
1347                                                 bset->dim, 0);
1348         bmap = add_constraints(ctx, bmap, bmap_domain, 0, 0);
1349
1350         bmap = isl_basic_map_simplify(ctx, bmap);
1351         return isl_basic_map_finalize(ctx, bmap);
1352 error:
1353         isl_basic_map_free(ctx, bmap);
1354         isl_basic_set_free(ctx, bset);
1355         return NULL;
1356 }
1357
1358 struct isl_basic_map *isl_basic_map_intersect_range(
1359                 struct isl_ctx *ctx, struct isl_basic_map *bmap,
1360                 struct isl_basic_set *bset)
1361 {
1362         struct isl_basic_map *bmap_range;
1363
1364         if (!bmap || !bset)
1365                 goto error;
1366
1367         isl_assert(ctx, bset->dim == bmap->n_out, goto error);
1368         isl_assert(ctx, bset->nparam == bmap->nparam, goto error);
1369
1370         bmap = isl_basic_map_extend(ctx, bmap,
1371                         bset->nparam, bmap->n_in, bmap->n_out,
1372                         bset->n_div, bset->n_eq, bset->n_ineq);
1373         if (!bmap)
1374                 goto error;
1375         bmap_range = isl_basic_map_from_basic_set(ctx, bset,
1376                                                 0, bset->dim);
1377         bmap = add_constraints(ctx, bmap, bmap_range, 0, 0);
1378
1379         bmap = isl_basic_map_simplify(ctx, bmap);
1380         return isl_basic_map_finalize(ctx, bmap);
1381 error:
1382         isl_basic_map_free(ctx, bmap);
1383         isl_basic_set_free(ctx, bset);
1384         return NULL;
1385 }
1386
1387 struct isl_basic_map *isl_basic_map_intersect(
1388                 struct isl_ctx *ctx, struct isl_basic_map *bmap1,
1389                 struct isl_basic_map *bmap2)
1390 {
1391         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
1392         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
1393         isl_assert(ctx, bmap1->n_out == bmap2->n_out, goto error);
1394
1395         bmap1 = isl_basic_map_extend(ctx, bmap1,
1396                         bmap1->nparam, bmap1->n_in, bmap1->n_out,
1397                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
1398         if (!bmap1)
1399                 goto error;
1400         bmap1 = add_constraints(ctx, bmap1, bmap2, 0, 0);
1401
1402         bmap1 = isl_basic_map_simplify(ctx, bmap1);
1403         return isl_basic_map_finalize(ctx, bmap1);
1404 error:
1405         isl_basic_map_free(ctx, bmap1);
1406         isl_basic_map_free(ctx, bmap2);
1407         return NULL;
1408 }
1409
1410 struct isl_basic_set *isl_basic_set_intersect(
1411                 struct isl_ctx *ctx, struct isl_basic_set *bset1,
1412                 struct isl_basic_set *bset2)
1413 {
1414         return (struct isl_basic_set *)
1415                 isl_basic_map_intersect(ctx,
1416                         (struct isl_basic_map *)bset1,
1417                         (struct isl_basic_map *)bset2);
1418 }
1419
1420 struct isl_map *isl_map_intersect(struct isl_ctx *ctx, struct isl_map *map1,
1421                 struct isl_map *map2)
1422 {
1423         unsigned flags = 0;
1424         struct isl_map *result;
1425         int i, j;
1426
1427         if (!map1 || !map2)
1428                 goto error;
1429
1430         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
1431             F_ISSET(map2, ISL_MAP_DISJOINT))
1432                 FL_SET(flags, ISL_MAP_DISJOINT);
1433
1434         result = isl_map_alloc(ctx, map1->nparam, map1->n_in, map1->n_out,
1435                                 map1->n * map2->n, flags);
1436         if (!result)
1437                 goto error;
1438         for (i = 0; i < map1->n; ++i)
1439                 for (j = 0; j < map2->n; ++j) {
1440                         struct isl_basic_map *part;
1441                         part = isl_basic_map_intersect(ctx,
1442                                     isl_basic_map_copy(ctx, map1->p[i]),
1443                                     isl_basic_map_copy(ctx, map2->p[j]));
1444                         if (isl_basic_map_is_empty(ctx, part))
1445                                 isl_basic_map_free(ctx, part);
1446                         else
1447                                 result = isl_map_add(ctx, result, part);
1448                         if (!result)
1449                                 goto error;
1450                 }
1451         isl_map_free(ctx, map1);
1452         isl_map_free(ctx, map2);
1453         return result;
1454 error:
1455         isl_map_free(ctx, map1);
1456         isl_map_free(ctx, map2);
1457         return NULL;
1458 }
1459
1460 struct isl_basic_map *isl_basic_map_reverse(struct isl_ctx *ctx,
1461                 struct isl_basic_map *bmap)
1462 {
1463         struct isl_basic_set *bset;
1464         unsigned in;
1465
1466         if (!bmap)
1467                 return NULL;
1468         bmap = isl_basic_map_cow(ctx, bmap);
1469         if (!bmap)
1470                 return NULL;
1471         in = bmap->n_in;
1472         bset = isl_basic_set_from_basic_map(ctx, bmap);
1473         bset = isl_basic_set_swap_vars(ctx, bset, in);
1474         return isl_basic_map_from_basic_set(ctx, bset, bset->dim-in, in);
1475 }
1476
1477 /* Turn final n dimensions into existentially quantified variables.
1478  */
1479 struct isl_basic_set *isl_basic_set_project_out(
1480                 struct isl_ctx *ctx, struct isl_basic_set *bset,
1481                 unsigned n, unsigned flags)
1482 {
1483         int i;
1484         size_t row_size;
1485         isl_int **new_div;
1486         isl_int *old;
1487
1488         if (!bset)
1489                 return NULL;
1490
1491         isl_assert(ctx, n <= bset->dim, goto error);
1492
1493         if (n == 0)
1494                 return bset;
1495
1496         bset = isl_basic_set_cow(ctx, bset);
1497
1498         row_size = 1 + bset->nparam + bset->dim + bset->extra;
1499         old = bset->block2.data;
1500         bset->block2 = isl_blk_extend(ctx, bset->block2,
1501                                         (bset->extra + n) * (1 + row_size));
1502         if (!bset->block2.data)
1503                 goto error;
1504         new_div = isl_alloc_array(ctx, isl_int *, bset->extra + n);
1505         if (!new_div)
1506                 goto error;
1507         for (i = 0; i < n; ++i) {
1508                 new_div[i] = bset->block2.data +
1509                                 (bset->extra + i) * (1 + row_size);
1510                 isl_seq_clr(new_div[i], 1 + row_size);
1511         }
1512         for (i = 0; i < bset->extra; ++i)
1513                 new_div[n + i] = bset->block2.data + (bset->div[i] - old);
1514         free(bset->div);
1515         bset->div = new_div;
1516         bset->n_div += n;
1517         bset->extra += n;
1518         bset->dim -= n;
1519         bset = isl_basic_set_simplify(ctx, bset);
1520         return isl_basic_set_finalize(ctx, bset);
1521 error:
1522         isl_basic_set_free(ctx, bset);
1523         return NULL;
1524 }
1525
1526 struct isl_basic_map *isl_basic_map_apply_range(
1527                 struct isl_ctx *ctx, struct isl_basic_map *bmap1,
1528                 struct isl_basic_map *bmap2)
1529 {
1530         struct isl_basic_set *bset;
1531         unsigned n_in, n_out;
1532
1533         if (!bmap1 || !bmap2)
1534                 goto error;
1535
1536         isl_assert(ctx, bmap1->n_out == bmap2->n_in, goto error);
1537         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
1538
1539         n_in = bmap1->n_in;
1540         n_out = bmap2->n_out;
1541
1542         bmap2 = isl_basic_map_reverse(ctx, bmap2);
1543         if (!bmap2)
1544                 goto error;
1545         bmap1 = isl_basic_map_extend(ctx, bmap1, bmap1->nparam,
1546                         bmap1->n_in + bmap2->n_in, bmap1->n_out,
1547                         bmap2->extra, bmap2->n_eq, bmap2->n_ineq);
1548         if (!bmap1)
1549                 goto error;
1550         bmap1 = add_constraints(ctx, bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
1551         bmap1 = isl_basic_map_simplify(ctx, bmap1);
1552         bset = isl_basic_set_from_basic_map(ctx, bmap1);
1553         bset = isl_basic_set_project_out(ctx, bset,
1554                                                 bset->dim - (n_in + n_out), 0);
1555         return isl_basic_map_from_basic_set(ctx, bset, n_in, n_out);
1556 error:
1557         isl_basic_map_free(ctx, bmap1);
1558         isl_basic_map_free(ctx, bmap2);
1559         return NULL;
1560 }
1561
1562 struct isl_basic_map *isl_basic_map_apply_domain(
1563                 struct isl_ctx *ctx, struct isl_basic_map *bmap1,
1564                 struct isl_basic_map *bmap2)
1565 {
1566         if (!bmap1 || !bmap2)
1567                 goto error;
1568
1569         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
1570         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
1571
1572         bmap1 = isl_basic_map_reverse(ctx, bmap1);
1573         bmap1 = isl_basic_map_apply_range(ctx, bmap1, bmap2);
1574         return isl_basic_map_reverse(ctx, bmap1);
1575 error:
1576         isl_basic_map_free(ctx, bmap1);
1577         isl_basic_map_free(ctx, bmap2);
1578         return NULL;
1579 }
1580
1581 static struct isl_basic_map *var_equal(struct isl_ctx *ctx,
1582                 struct isl_basic_map *bmap, unsigned pos)
1583 {
1584         int i;
1585         i = isl_basic_map_alloc_equality(ctx, bmap);
1586         if (i < 0)
1587                 goto error;
1588         isl_seq_clr(bmap->eq[i],
1589                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1590         isl_int_set_si(bmap->eq[i][1+bmap->nparam+pos], -1);
1591         isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1592         return isl_basic_map_finalize(ctx, bmap);
1593 error:
1594         isl_basic_map_free(ctx, bmap);
1595         return NULL;
1596 }
1597
1598 static struct isl_basic_map *var_more(struct isl_ctx *ctx,
1599                 struct isl_basic_map *bmap, unsigned pos)
1600 {
1601         int i;
1602         i = isl_basic_map_alloc_inequality(ctx, bmap);
1603         if (i < 0)
1604                 goto error;
1605         isl_seq_clr(bmap->ineq[i],
1606                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1607         isl_int_set_si(bmap->ineq[i][0], -1);
1608         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], -1);
1609         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1610         return isl_basic_map_finalize(ctx, bmap);
1611 error:
1612         isl_basic_map_free(ctx, bmap);
1613         return NULL;
1614 }
1615
1616 static struct isl_basic_map *var_less(struct isl_ctx *ctx,
1617                 struct isl_basic_map *bmap, unsigned pos)
1618 {
1619         int i;
1620         i = isl_basic_map_alloc_inequality(ctx, bmap);
1621         if (i < 0)
1622                 goto error;
1623         isl_seq_clr(bmap->ineq[i],
1624                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1625         isl_int_set_si(bmap->ineq[i][0], -1);
1626         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], 1);
1627         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], -1);
1628         return isl_basic_map_finalize(ctx, bmap);
1629 error:
1630         isl_basic_map_free(ctx, bmap);
1631         return NULL;
1632 }
1633
1634 struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
1635                 unsigned nparam, unsigned in, unsigned out, unsigned n_equal)
1636 {
1637         int i;
1638         struct isl_basic_map *bmap;
1639         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, n_equal, 0);
1640         if (!bmap)
1641                 return NULL;
1642         for (i = 0; i < n_equal && bmap; ++i)
1643                 bmap = var_equal(ctx, bmap, i);
1644         return isl_basic_map_finalize(ctx, bmap);
1645 }
1646
1647 struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
1648                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1649 {
1650         int i;
1651         struct isl_basic_map *bmap;
1652         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1653         if (!bmap)
1654                 return NULL;
1655         for (i = 0; i < pos && bmap; ++i)
1656                 bmap = var_equal(ctx, bmap, i);
1657         if (bmap)
1658                 bmap = var_less(ctx, bmap, pos);
1659         return isl_basic_map_finalize(ctx, bmap);
1660 }
1661
1662 struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
1663                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1664 {
1665         int i;
1666         struct isl_basic_map *bmap;
1667         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1668         if (!bmap)
1669                 return NULL;
1670         for (i = 0; i < pos && bmap; ++i)
1671                 bmap = var_equal(ctx, bmap, i);
1672         if (bmap)
1673                 bmap = var_more(ctx, bmap, pos);
1674         return isl_basic_map_finalize(ctx, bmap);
1675 }
1676
1677 struct isl_basic_map *isl_basic_map_from_basic_set(
1678                 struct isl_ctx *ctx, struct isl_basic_set *bset,
1679                 unsigned n_in, unsigned n_out)
1680 {
1681         struct isl_basic_map *bmap;
1682
1683         bset = isl_basic_set_cow(ctx, bset);
1684         if (!bset)
1685                 return NULL;
1686
1687         isl_assert(ctx, bset->dim == n_in + n_out, goto error);
1688         bmap = (struct isl_basic_map *) bset;
1689         bmap->n_in = n_in;
1690         bmap->n_out = n_out;
1691         return isl_basic_map_finalize(ctx, bmap);
1692 error:
1693         isl_basic_set_free(ctx, bset);
1694         return NULL;
1695 }
1696
1697 struct isl_basic_set *isl_basic_set_from_basic_map(
1698                 struct isl_ctx *ctx, struct isl_basic_map *bmap)
1699 {
1700         if (!bmap)
1701                 goto error;
1702         if (bmap->n_in == 0)
1703                 return (struct isl_basic_set *)bmap;
1704         bmap = isl_basic_map_cow(ctx, bmap);
1705         if (!bmap)
1706                 goto error;
1707         bmap->n_out += bmap->n_in;
1708         bmap->n_in = 0;
1709         bmap = isl_basic_map_finalize(ctx, bmap);
1710         return (struct isl_basic_set *)bmap;
1711 error:
1712         return NULL;
1713 }
1714
1715 static struct isl_basic_set *isl_basic_map_underlying_set(
1716                 struct isl_ctx *ctx, struct isl_basic_map *bmap)
1717 {
1718         if (!bmap)
1719                 goto error;
1720         if (bmap->nparam == 0 && bmap->n_in == 0 && bmap->n_div == 0)
1721                 return (struct isl_basic_set *)bmap;
1722         bmap = isl_basic_map_cow(ctx, bmap);
1723         if (!bmap)
1724                 goto error;
1725         bmap->n_out += bmap->nparam + bmap->n_in + bmap->n_div;
1726         bmap->nparam = 0;
1727         bmap->n_in = 0;
1728         bmap->n_div = 0;
1729         bmap = isl_basic_map_finalize(ctx, bmap);
1730         return (struct isl_basic_set *)bmap;
1731 error:
1732         return NULL;
1733 }
1734
1735 struct isl_basic_set *isl_basic_map_domain(struct isl_ctx *ctx,
1736                 struct isl_basic_map *bmap)
1737 {
1738         struct isl_basic_set *domain;
1739         unsigned n_out;
1740         if (!bmap)
1741                 return NULL;
1742         n_out = bmap->n_out;
1743         domain = isl_basic_set_from_basic_map(ctx, bmap);
1744         return isl_basic_set_project_out(ctx, domain, n_out, 0);
1745 }
1746
1747 struct isl_basic_set *isl_basic_map_range(struct isl_ctx *ctx,
1748                 struct isl_basic_map *bmap)
1749 {
1750         return isl_basic_map_domain(ctx,
1751                         isl_basic_map_reverse(ctx, bmap));
1752 }
1753
1754 struct isl_set *isl_map_range(struct isl_ctx *ctx, struct isl_map *map)
1755 {
1756         int i;
1757         struct isl_set *set;
1758
1759         if (!map)
1760                 goto error;
1761         map = isl_map_cow(ctx, map);
1762         if (!map)
1763                 goto error;
1764
1765         set = (struct isl_set *) map;
1766         set->zero = 0;
1767         for (i = 0; i < map->n; ++i) {
1768                 set->p[i] = isl_basic_map_range(ctx, map->p[i]);
1769                 if (!set->p[i])
1770                         goto error;
1771         }
1772         F_CLR(set, ISL_MAP_DISJOINT);
1773         return set;
1774 error:
1775         isl_map_free(ctx, map);
1776         return NULL;
1777 }
1778
1779 struct isl_map *isl_map_from_set(
1780                 struct isl_ctx *ctx, struct isl_set *set,
1781                 unsigned n_in, unsigned n_out)
1782 {
1783         int i;
1784         struct isl_map *map = NULL;
1785
1786         if (!set)
1787                 return NULL;
1788         isl_assert(ctx, set->dim == n_in + n_out, goto error);
1789         set = isl_set_cow(ctx, set);
1790         if (!set)
1791                 return NULL;
1792         map = (struct isl_map *)set;
1793         for (i = 0; i < set->n; ++i) {
1794                 map->p[i] = isl_basic_map_from_basic_set(ctx,
1795                                 set->p[i], n_in, n_out);
1796                 if (!map->p[i])
1797                         goto error;
1798         }
1799         map->n_in = n_in;
1800         map->n_out = n_out;
1801         return map;
1802 error:
1803         isl_set_free(ctx, set);
1804         return NULL;
1805 }
1806
1807 struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
1808                 unsigned nparam, unsigned in, unsigned out, int n,
1809                 unsigned flags)
1810 {
1811         struct isl_map *map;
1812
1813         map = isl_alloc(ctx, struct isl_map,
1814                         sizeof(struct isl_map) +
1815                         n * sizeof(struct isl_basic_map *));
1816         if (!map)
1817                 return NULL;
1818
1819         map->ref = 1;
1820         map->size = n;
1821         map->n = 0;
1822         map->nparam = nparam;
1823         map->n_in = in;
1824         map->n_out = out;
1825         map->flags = flags;
1826         return map;
1827 }
1828
1829 struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
1830                 unsigned nparam, unsigned in, unsigned out)
1831 {
1832         struct isl_basic_map *bmap;
1833         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 1, 0);
1834         bmap = isl_basic_map_set_to_empty(ctx, bmap);
1835         return bmap;
1836 }
1837
1838 struct isl_map *isl_map_empty(struct isl_ctx *ctx,
1839                 unsigned nparam, unsigned in, unsigned out)
1840 {
1841         return isl_map_alloc(ctx, nparam, in, out, 0, ISL_MAP_DISJOINT);
1842 }
1843
1844 struct isl_set *isl_set_empty(struct isl_ctx *ctx,
1845                 unsigned nparam, unsigned dim)
1846 {
1847         return isl_set_alloc(ctx, nparam, dim, 0, ISL_MAP_DISJOINT);
1848 }
1849
1850 struct isl_map *isl_map_dup(struct isl_ctx *ctx, struct isl_map *map)
1851 {
1852         int i;
1853         struct isl_map *dup;
1854
1855         dup = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out, map->n,
1856                                 map->flags);
1857         for (i = 0; i < map->n; ++i)
1858                 dup = isl_map_add(ctx, dup,
1859                                   isl_basic_map_dup(ctx, map->p[i]));
1860         return dup;
1861 }
1862
1863 struct isl_map *isl_map_add(struct isl_ctx *ctx, struct isl_map *map,
1864                                 struct isl_basic_map *bmap)
1865 {
1866         if (!bmap || !map)
1867                 goto error;
1868         isl_assert(ctx, map->nparam == bmap->nparam, goto error);
1869         isl_assert(ctx, map->n_in == bmap->n_in, goto error);
1870         isl_assert(ctx, map->n_out == bmap->n_out, goto error);
1871         isl_assert(ctx, map->n < map->size, goto error);
1872         map->p[map->n] = bmap;
1873         map->n++;
1874         return map;
1875 error:
1876         if (map)
1877                 isl_map_free(ctx, map);
1878         if (bmap)
1879                 isl_basic_map_free(ctx, bmap);
1880         return NULL;
1881 }
1882
1883 void isl_map_free(struct isl_ctx *ctx, struct isl_map *map)
1884 {
1885         int i;
1886
1887         if (!map)
1888                 return;
1889
1890         if (--map->ref > 0)
1891                 return;
1892
1893         for (i = 0; i < map->n; ++i)
1894                 isl_basic_map_free(ctx, map->p[i]);
1895         free(map);
1896 }
1897
1898 struct isl_map *isl_map_extend(struct isl_ctx *ctx, struct isl_map *base,
1899                 unsigned nparam, unsigned n_in, unsigned n_out)
1900 {
1901         int i;
1902
1903         base = isl_map_cow(ctx, base);
1904         if (!base)
1905                 return NULL;
1906
1907         isl_assert(ctx, base->nparam <= nparam, goto error);
1908         isl_assert(ctx, base->n_in <= n_in, goto error);
1909         isl_assert(ctx, base->n_out <= n_out, goto error);
1910         base->nparam = nparam;
1911         base->n_in = n_in;
1912         base->n_out = n_out;
1913         for (i = 0; i < base->n; ++i) {
1914                 base->p[i] = isl_basic_map_extend(ctx, base->p[i],
1915                                 nparam, n_in, n_out, 0, 0, 0);
1916                 if (!base->p[i])
1917                         goto error;
1918         }
1919         return base;
1920 error:
1921         isl_map_free(ctx, base);
1922         return NULL;
1923 }
1924
1925 struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_ctx *ctx,
1926                 struct isl_basic_map *bmap,
1927                 unsigned input, int value)
1928 {
1929         int j;
1930
1931         bmap = isl_basic_map_cow(ctx, bmap);
1932         if (!bmap)
1933                 return NULL;
1934         isl_assert(ctx, input < bmap->n_in, goto error);
1935
1936         bmap = isl_basic_map_extend(ctx, bmap,
1937                         bmap->nparam, bmap->n_in, bmap->n_out, 0, 1, 0);
1938         j = isl_basic_map_alloc_equality(ctx, bmap);
1939         if (j < 0)
1940                 goto error;
1941         isl_int_set_si(bmap->eq[j][1+bmap->nparam+input], -1);
1942         isl_int_set_si(bmap->eq[j][0], value);
1943         bmap = isl_basic_map_simplify(ctx, bmap);
1944         return isl_basic_map_finalize(ctx, bmap);
1945 error:
1946         isl_basic_map_free(ctx, bmap);
1947         return NULL;
1948 }
1949
1950 struct isl_map *isl_map_fix_input_si(struct isl_ctx *ctx, struct isl_map *map,
1951                 unsigned input, int value)
1952 {
1953         int i;
1954
1955         map = isl_map_cow(ctx, map);
1956         if (!map)
1957                 return NULL;
1958
1959         isl_assert(ctx, input < map->n_in, goto error);
1960         for (i = 0; i < map->n; ++i) {
1961                 map->p[i] = isl_basic_map_fix_input_si(ctx, map->p[i],
1962                                                                 input, value);
1963                 if (!map->p[i])
1964                         goto error;
1965         }
1966         return map;
1967 error:
1968         isl_map_free(ctx, map);
1969         return NULL;
1970 }
1971
1972 struct isl_map *isl_map_reverse(struct isl_ctx *ctx, struct isl_map *map)
1973 {
1974         int i;
1975         unsigned t;
1976
1977         map = isl_map_cow(ctx, map);
1978         if (!map)
1979                 return NULL;
1980
1981         t = map->n_in;
1982         map->n_in = map->n_out;
1983         map->n_out = t;
1984         for (i = 0; i < map->n; ++i) {
1985                 map->p[i] = isl_basic_map_reverse(ctx, map->p[i]);
1986                 if (!map->p[i])
1987                         goto error;
1988         }
1989         return map;
1990 error:
1991         isl_map_free(ctx, map);
1992         return NULL;
1993 }
1994
1995 struct isl_map *isl_basic_map_lexmax(struct isl_ctx *ctx,
1996                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
1997                 struct isl_set **empty)
1998 {
1999         return isl_pip_basic_map_lexmax(ctx, bmap, dom, empty);
2000 }
2001
2002 struct isl_map *isl_basic_map_lexmin(struct isl_ctx *ctx,
2003                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2004                 struct isl_set **empty)
2005 {
2006         return isl_pip_basic_map_lexmin(ctx, bmap, dom, empty);
2007 }
2008
2009 struct isl_set *isl_basic_set_lexmin(struct isl_ctx *ctx,
2010                 struct isl_basic_set *bset)
2011 {
2012         struct isl_basic_map *bmap = NULL;
2013         struct isl_basic_set *dom = NULL;
2014         struct isl_map *min;
2015
2016         if (!bset)
2017                 goto error;
2018         bmap = isl_basic_map_from_basic_set(ctx, bset, 0, bset->dim);
2019         if (!bmap)
2020                 goto error;
2021         dom = isl_basic_set_alloc(ctx, bmap->nparam, 0, 0, 0, 0);
2022         if (!dom)
2023                 goto error;
2024         min = isl_basic_map_lexmin(ctx, bmap, dom, NULL);
2025         return isl_map_range(ctx, min);
2026 error:
2027         isl_basic_map_free(ctx, bmap);
2028         return NULL;
2029 }
2030
2031 struct isl_map *isl_basic_map_compute_divs(struct isl_ctx *ctx,
2032                 struct isl_basic_map *bmap)
2033 {
2034         if (bmap->n_div == 0)
2035                 return isl_map_from_basic_map(ctx, bmap);
2036         return isl_pip_basic_map_compute_divs(ctx, bmap);
2037 }
2038
2039 struct isl_map *isl_map_compute_divs(struct isl_ctx *ctx, struct isl_map *map)
2040 {
2041         int i;
2042         struct isl_map *res;
2043
2044         if (!map)
2045                 return NULL;
2046         if (map->n == 0)
2047                 return map;
2048         res = isl_basic_map_compute_divs(ctx,
2049                 isl_basic_map_copy(ctx, map->p[0]));
2050         for (i = 1 ; i < map->n; ++i) {
2051                 struct isl_map *r2;
2052                 r2 = isl_basic_map_compute_divs(ctx,
2053                         isl_basic_map_copy(ctx, map->p[i]));
2054                 if (F_ISSET(map, ISL_MAP_DISJOINT))
2055                         res = isl_map_union_disjoint(ctx, res, r2);
2056                 else
2057                         res = isl_map_union(ctx, res, r2);
2058         }
2059         isl_map_free(ctx, map);
2060
2061         return res;
2062 }
2063
2064 struct isl_set *isl_basic_set_compute_divs(struct isl_ctx *ctx,
2065                 struct isl_basic_set *bset)
2066 {
2067         return (struct isl_set *)
2068                 isl_basic_map_compute_divs(ctx,
2069                         (struct isl_basic_map *)bset);
2070 }
2071
2072 struct isl_set *isl_map_domain(struct isl_ctx *ctx, struct isl_map *map)
2073 {
2074         int i;
2075         struct isl_set *set;
2076
2077         if (!map)
2078                 goto error;
2079
2080         map = isl_map_cow(ctx, map);
2081         if (!map)
2082                 return NULL;
2083
2084         set = (struct isl_set *)map;
2085         set->dim = map->n_in;
2086         set->zero = 0;
2087         for (i = 0; i < map->n; ++i) {
2088                 set->p[i] = isl_basic_map_domain(ctx, map->p[i]);
2089                 if (!set->p[i])
2090                         goto error;
2091         }
2092         F_CLR(set, ISL_MAP_DISJOINT);
2093         return set;
2094 error:
2095         isl_map_free(ctx, map);
2096         return NULL;
2097 }
2098
2099 struct isl_map *isl_map_union_disjoint(struct isl_ctx *ctx,
2100                         struct isl_map *map1, struct isl_map *map2)
2101 {
2102         int i;
2103         unsigned flags = 0;
2104         struct isl_map *map = NULL;
2105
2106         if (!map1 || !map2)
2107                 goto error;
2108
2109         if (map1->n == 0) {
2110                 isl_map_free(ctx, map1);
2111                 return map2;
2112         }
2113         if (map2->n == 0) {
2114                 isl_map_free(ctx, map2);
2115                 return map1;
2116         }
2117
2118         isl_assert(ctx, map1->nparam == map2->nparam, goto error);
2119         isl_assert(ctx, map1->n_in == map2->n_in, goto error);
2120         isl_assert(ctx, map1->n_out == map2->n_out, goto error);
2121
2122         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
2123             F_ISSET(map2, ISL_MAP_DISJOINT))
2124                 FL_SET(flags, ISL_MAP_DISJOINT);
2125
2126         map = isl_map_alloc(ctx, map1->nparam, map1->n_in, map1->n_out,
2127                                 map1->n + map2->n, flags);
2128         if (!map)
2129                 goto error;
2130         for (i = 0; i < map1->n; ++i) {
2131                 map = isl_map_add(ctx, map,
2132                                   isl_basic_map_copy(ctx, map1->p[i]));
2133                 if (!map)
2134                         goto error;
2135         }
2136         for (i = 0; i < map2->n; ++i) {
2137                 map = isl_map_add(ctx, map,
2138                                   isl_basic_map_copy(ctx, map2->p[i]));
2139                 if (!map)
2140                         goto error;
2141         }
2142         isl_map_free(ctx, map1);
2143         isl_map_free(ctx, map2);
2144         return map;
2145 error:
2146         isl_map_free(ctx, map);
2147         isl_map_free(ctx, map1);
2148         isl_map_free(ctx, map2);
2149         return NULL;
2150 }
2151
2152 struct isl_map *isl_map_union(struct isl_ctx *ctx,
2153                         struct isl_map *map1, struct isl_map *map2)
2154 {
2155         map1 = isl_map_union_disjoint(ctx, map1, map2);
2156         if (!map1)
2157                 return NULL;
2158         if (map1->n > 1)
2159                 F_CLR(map1, ISL_MAP_DISJOINT);
2160         return map1;
2161 }
2162
2163 struct isl_set *isl_set_union_disjoint(struct isl_ctx *ctx,
2164                         struct isl_set *set1, struct isl_set *set2)
2165 {
2166         return (struct isl_set *)
2167                 isl_map_union_disjoint(ctx,
2168                         (struct isl_map *)set1, (struct isl_map *)set2);
2169 }
2170
2171 struct isl_set *isl_set_union(struct isl_ctx *ctx,
2172                         struct isl_set *set1, struct isl_set *set2)
2173 {
2174         return (struct isl_set *)
2175                 isl_map_union(ctx,
2176                         (struct isl_map *)set1, (struct isl_map *)set2);
2177 }
2178
2179 struct isl_map *isl_map_intersect_range(
2180                 struct isl_ctx *ctx, struct isl_map *map,
2181                 struct isl_set *set)
2182 {
2183         unsigned flags = 0;
2184         struct isl_map *result;
2185         int i, j;
2186
2187         if (!map || !set)
2188                 goto error;
2189
2190         if (F_ISSET(map, ISL_MAP_DISJOINT) &&
2191             F_ISSET(set, ISL_MAP_DISJOINT))
2192                 FL_SET(flags, ISL_MAP_DISJOINT);
2193
2194         result = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out,
2195                                 map->n * set->n, flags);
2196         if (!result)
2197                 goto error;
2198         for (i = 0; i < map->n; ++i)
2199                 for (j = 0; j < set->n; ++j) {
2200                         result = isl_map_add(ctx, result,
2201                             isl_basic_map_intersect_range(ctx,
2202                                 isl_basic_map_copy(ctx, map->p[i]),
2203                                 isl_basic_set_copy(ctx, set->p[j])));
2204                         if (!result)
2205                                 goto error;
2206                 }
2207         isl_map_free(ctx, map);
2208         isl_set_free(ctx, set);
2209         return result;
2210 error:
2211         isl_map_free(ctx, map);
2212         isl_set_free(ctx, set);
2213         return NULL;
2214 }
2215
2216 struct isl_map *isl_map_intersect_domain(
2217                 struct isl_ctx *ctx, struct isl_map *map,
2218                 struct isl_set *set)
2219 {
2220         return isl_map_reverse(ctx,
2221                 isl_map_intersect_range(ctx, isl_map_reverse(ctx, map), set));
2222 }
2223
2224 struct isl_map *isl_map_apply_domain(
2225                 struct isl_ctx *ctx, struct isl_map *map1,
2226                 struct isl_map *map2)
2227 {
2228         if (!map1 || !map2)
2229                 goto error;
2230         map1 = isl_map_reverse(ctx, map1);
2231         map1 = isl_map_apply_range(ctx, map1, map2);
2232         return isl_map_reverse(ctx, map1);
2233 error:
2234         isl_map_free(ctx, map1);
2235         isl_map_free(ctx, map2);
2236         return NULL;
2237 }
2238
2239 struct isl_map *isl_map_apply_range(
2240                 struct isl_ctx *ctx, struct isl_map *map1,
2241                 struct isl_map *map2)
2242 {
2243         struct isl_map *result;
2244         int i, j;
2245
2246         if (!map1 || !map2)
2247                 goto error;
2248
2249         isl_assert(ctx, map1->nparam == map2->nparam, goto error);
2250         isl_assert(ctx, map1->n_out == map2->n_in, goto error);
2251
2252         result = isl_map_alloc(ctx, map1->nparam, map1->n_in, map2->n_out,
2253                                 map1->n * map2->n, 0);
2254         if (!result)
2255                 goto error;
2256         for (i = 0; i < map1->n; ++i)
2257                 for (j = 0; j < map2->n; ++j) {
2258                         result = isl_map_add(ctx, result,
2259                             isl_basic_map_apply_range(ctx,
2260                                 isl_basic_map_copy(ctx, map1->p[i]),
2261                                 isl_basic_map_copy(ctx, map2->p[j])));
2262                         if (!result)
2263                                 goto error;
2264                 }
2265         isl_map_free(ctx, map1);
2266         isl_map_free(ctx, map2);
2267         return result;
2268 error:
2269         isl_map_free(ctx, map1);
2270         isl_map_free(ctx, map2);
2271         return NULL;
2272 }
2273
2274 /*
2275  * returns range - domain
2276  */
2277 struct isl_basic_set *isl_basic_map_deltas(struct isl_ctx *ctx,
2278                                 struct isl_basic_map *bmap)
2279 {
2280         struct isl_basic_set *bset;
2281         unsigned dim;
2282         int i;
2283
2284         if (!bmap)
2285                 goto error;
2286         isl_assert(ctx, bmap->n_in == bmap->n_out, goto error);
2287         dim = bmap->n_in;
2288         bset = isl_basic_set_from_basic_map(ctx, bmap);
2289         bset = isl_basic_set_extend(ctx, bset, bset->nparam, 3*dim, 0,
2290                                         dim, 0);
2291         bset = isl_basic_set_swap_vars(ctx, bset, 2*dim);
2292         for (i = 0; i < dim; ++i) {
2293                 int j = isl_basic_map_alloc_equality(ctx,
2294                                             (struct isl_basic_map *)bset);
2295                 if (j < 0)
2296                         goto error;
2297                 isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
2298                 isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
2299                 isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
2300         }
2301         return isl_basic_set_project_out(ctx, bset, 2*dim, 0);
2302 error:
2303         isl_basic_map_free(ctx, bmap);
2304         return NULL;
2305 }
2306
2307 /* If the only constraints a div d=floor(f/m)
2308  * appears in are its two defining constraints
2309  *
2310  *      f - m d >=0
2311  *      -(f - (m - 1)) + m d >= 0
2312  *
2313  * then it can safely be removed.
2314  */
2315 static int div_is_redundant(struct isl_ctx *ctx, struct isl_basic_map *bmap,
2316                                 int div)
2317 {
2318         int i;
2319         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2320
2321         for (i = 0; i < bmap->n_eq; ++i)
2322                 if (!isl_int_is_zero(bmap->eq[i][pos]))
2323                         return 0;
2324
2325         for (i = 0; i < bmap->n_ineq; ++i) {
2326                 if (isl_int_is_zero(bmap->ineq[i][pos]))
2327                         continue;
2328                 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2329                         int neg;
2330                         isl_int_sub(bmap->div[div][1],
2331                                         bmap->div[div][1], bmap->div[div][0]);
2332                         isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
2333                         neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
2334                         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
2335                         isl_int_add(bmap->div[div][1],
2336                                         bmap->div[div][1], bmap->div[div][0]);
2337                         if (!neg)
2338                                 return 0;
2339                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2340                                                     bmap->n_div-div-1) != -1)
2341                                 return 0;
2342                 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2343                         if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
2344                                 return 0;
2345                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2346                                                     bmap->n_div-div-1) != -1)
2347                                 return 0;
2348                 } else
2349                         return 0;
2350         }
2351
2352         for (i = 0; i < bmap->n_div; ++i)
2353                 if (!isl_int_is_zero(bmap->div[i][1+pos]))
2354                         return 0;
2355
2356         return 1;
2357 }
2358
2359 /*
2360  * Remove divs that don't occur in any of the constraints or other divs.
2361  * These can arise when dropping some of the variables in a quast
2362  * returned by piplib.
2363  */
2364 static struct isl_basic_map *remove_redundant_divs(struct isl_ctx *ctx,
2365                 struct isl_basic_map *bmap)
2366 {
2367         int i;
2368
2369         for (i = bmap->n_div-1; i >= 0; --i) {
2370                 if (!div_is_redundant(ctx, bmap, i))
2371                         continue;
2372                 bmap = isl_basic_map_drop_div(ctx, bmap, i);
2373         }
2374         return bmap;
2375 }
2376
2377 struct isl_basic_map *isl_basic_map_finalize(struct isl_ctx *ctx,
2378                 struct isl_basic_map *bmap)
2379 {
2380         bmap = remove_redundant_divs(ctx, bmap);
2381         if (!bmap)
2382                 return NULL;
2383         F_SET(bmap, ISL_BASIC_SET_FINAL);
2384         return bmap;
2385 }
2386
2387 struct isl_basic_set *isl_basic_set_finalize(struct isl_ctx *ctx,
2388                 struct isl_basic_set *bset)
2389 {
2390         return (struct isl_basic_set *)
2391                 isl_basic_map_finalize(ctx,
2392                         (struct isl_basic_map *)bset);
2393 }
2394
2395 struct isl_set *isl_set_finalize(struct isl_ctx *ctx, struct isl_set *set)
2396 {
2397         int i;
2398
2399         if (!set)
2400                 return NULL;
2401         for (i = 0; i < set->n; ++i) {
2402                 set->p[i] = isl_basic_set_finalize(ctx, set->p[i]);
2403                 if (!set->p[i])
2404                         goto error;
2405         }
2406         return set;
2407 error:
2408         isl_set_free(ctx, set);
2409         return NULL;
2410 }
2411
2412 struct isl_map *isl_map_finalize(struct isl_ctx *ctx, struct isl_map *map)
2413 {
2414         int i;
2415
2416         if (!map)
2417                 return NULL;
2418         for (i = 0; i < map->n; ++i) {
2419                 map->p[i] = isl_basic_map_finalize(ctx, map->p[i]);
2420                 if (!map->p[i])
2421                         goto error;
2422         }
2423         return map;
2424 error:
2425         isl_map_free(ctx, map);
2426         return NULL;
2427 }
2428
2429 struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
2430                 unsigned nparam, unsigned dim)
2431 {
2432         struct isl_basic_map *bmap;
2433         int i;
2434
2435         bmap = isl_basic_map_alloc(ctx, nparam, dim, dim, 0, dim, 0);
2436         if (!bmap)
2437                 goto error;
2438
2439         for (i = 0; i < dim; ++i) {
2440                 int j = isl_basic_map_alloc_equality(ctx, bmap);
2441                 if (j < 0)
2442                         goto error;
2443                 isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
2444                 isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
2445         }
2446         return isl_basic_map_finalize(ctx, bmap);
2447 error:
2448         isl_basic_map_free(ctx, bmap);
2449         return NULL;
2450 }
2451
2452 struct isl_map *isl_map_identity(struct isl_ctx *ctx,
2453                 unsigned nparam, unsigned dim)
2454 {
2455         struct isl_map *map = isl_map_alloc(ctx, nparam, dim, dim, 1,
2456                                                 ISL_MAP_DISJOINT);
2457         if (!map)
2458                 goto error;
2459         map = isl_map_add(ctx, map,
2460                 isl_basic_map_identity(ctx, nparam, dim));
2461         return map;
2462 error:
2463         isl_map_free(ctx, map);
2464         return NULL;
2465 }
2466
2467 int isl_set_is_equal(struct isl_ctx *ctx,
2468                 struct isl_set *set1, struct isl_set *set2)
2469 {
2470         return isl_map_is_equal(ctx,
2471                         (struct isl_map *)set1, (struct isl_map *)set2);
2472 }
2473
2474 int isl_set_is_subset(struct isl_ctx *ctx,
2475                 struct isl_set *set1, struct isl_set *set2)
2476 {
2477         return isl_map_is_subset(ctx,
2478                         (struct isl_map *)set1, (struct isl_map *)set2);
2479 }
2480
2481 int isl_basic_map_is_subset(struct isl_ctx *ctx,
2482                 struct isl_basic_map *bmap1,
2483                 struct isl_basic_map *bmap2)
2484 {
2485         int is_subset;
2486         struct isl_map *map1;
2487         struct isl_map *map2;
2488
2489         if (!bmap1 || !bmap2)
2490                 return -1;
2491
2492         map1 = isl_map_from_basic_map(ctx,
2493                         isl_basic_map_copy(ctx, bmap1));
2494         map2 = isl_map_from_basic_map(ctx,
2495                         isl_basic_map_copy(ctx, bmap2));
2496
2497         is_subset = isl_map_is_subset(ctx, map1, map2);
2498
2499         isl_map_free(ctx, map1);
2500         isl_map_free(ctx, map2);
2501
2502         return is_subset;
2503 }
2504
2505 int isl_map_is_empty(struct isl_ctx *ctx, struct isl_map *map)
2506 {
2507         int i;
2508         int is_empty;
2509
2510         if (!map)
2511                 return -1;
2512         for (i = 0; i < map->n; ++i) {
2513                 is_empty = isl_basic_map_is_empty(ctx, map->p[i]);
2514                 if (is_empty < 0)
2515                         return -1;
2516                 if (!is_empty)
2517                         return 0;
2518         }
2519         return 1;
2520 }
2521
2522 int isl_set_is_empty(struct isl_ctx *ctx, struct isl_set *set)
2523 {
2524         return isl_map_is_empty(ctx, (struct isl_map *)set);
2525 }
2526
2527 int isl_map_is_subset(struct isl_ctx *ctx, struct isl_map *map1,
2528                 struct isl_map *map2)
2529 {
2530         int i;
2531         int is_subset = 0;
2532         struct isl_map *diff;
2533
2534         if (!map1 || !map2)
2535                 return -1;
2536
2537         if (isl_map_is_empty(ctx, map1))
2538                 return 1;
2539
2540         if (isl_map_is_empty(ctx, map2))
2541                 return 0;
2542
2543         diff = isl_map_subtract(ctx, isl_map_copy(ctx, map1),
2544                                      isl_map_copy(ctx, map2));
2545         if (!diff)
2546                 return -1;
2547
2548         is_subset = isl_map_is_empty(ctx, diff);
2549         isl_map_free(ctx, diff);
2550
2551         return is_subset;
2552 }
2553
2554 int isl_map_is_equal(struct isl_ctx *ctx,
2555                 struct isl_map *map1, struct isl_map *map2)
2556 {
2557         int is_subset;
2558
2559         if (!map1 || !map2)
2560                 return -1;
2561         is_subset = isl_map_is_subset(ctx, map1, map2);
2562         if (is_subset != 1)
2563                 return is_subset;
2564         is_subset = isl_map_is_subset(ctx, map2, map1);
2565         return is_subset;
2566 }
2567
2568 int isl_basic_map_is_strict_subset(struct isl_ctx *ctx,
2569                 struct isl_basic_map *bmap1,
2570                 struct isl_basic_map *bmap2)
2571 {
2572         int is_subset;
2573
2574         if (!bmap1 || !bmap2)
2575                 return -1;
2576         is_subset = isl_basic_map_is_subset(ctx, bmap1, bmap2);
2577         if (is_subset != 1)
2578                 return is_subset;
2579         is_subset = isl_basic_map_is_subset(ctx, bmap2, bmap1);
2580         if (is_subset == -1)
2581                 return is_subset;
2582         return !is_subset;
2583 }
2584
2585 int isl_basic_map_is_empty(struct isl_ctx *ctx,
2586                 struct isl_basic_map *bmap)
2587 {
2588         struct isl_basic_set *bset = NULL;
2589         struct isl_vec *sample = NULL;
2590         int empty;
2591
2592         if (!bmap)
2593                 return -1;
2594
2595         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
2596                 return 1;
2597
2598         bset = isl_basic_map_underlying_set(ctx,
2599                         isl_basic_map_copy(ctx, bmap));
2600         if (!bset)
2601                 return -1;
2602         sample = isl_basic_set_sample(ctx, bset);
2603         if (!sample)
2604                 return -1;
2605         empty = sample->size == 0;
2606         isl_vec_free(ctx, sample);
2607
2608         return empty;
2609 }
2610
2611 struct isl_map *isl_basic_map_union(
2612                 struct isl_ctx *ctx, struct isl_basic_map *bmap1,
2613                 struct isl_basic_map *bmap2)
2614 {
2615         struct isl_map *map;
2616         if (!bmap1 || !bmap2)
2617                 return NULL;
2618
2619         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
2620         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
2621         isl_assert(ctx, bmap1->n_out == bmap2->n_out, goto error);
2622
2623         map = isl_map_alloc(ctx, bmap1->nparam,
2624                                 bmap1->n_in, bmap1->n_out, 2, 0);
2625         if (!map)
2626                 goto error;
2627         map = isl_map_add(ctx, map, bmap1);
2628         map = isl_map_add(ctx, map, bmap2);
2629         return map;
2630 error:
2631         isl_basic_map_free(ctx, bmap1);
2632         isl_basic_map_free(ctx, bmap2);
2633         return NULL;
2634 }
2635
2636 struct isl_set *isl_basic_set_union(
2637                 struct isl_ctx *ctx, struct isl_basic_set *bset1,
2638                 struct isl_basic_set *bset2)
2639 {
2640         return (struct isl_set *)isl_basic_map_union(ctx,
2641                                             (struct isl_basic_map *)bset1,
2642                                             (struct isl_basic_map *)bset2);
2643 }
2644
2645 /* Order divs such that any div only depends on previous divs */
2646 static struct isl_basic_map *order_divs(struct isl_ctx *ctx,
2647                 struct isl_basic_map *bmap)
2648 {
2649         int i;
2650         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
2651
2652         for (i = 0; i < bmap->n_div; ++i) {
2653                 int pos;
2654                 pos = isl_seq_first_non_zero(bmap->div[i]+1+1+off+i,
2655                                                             bmap->n_div-i);
2656                 if (pos == -1)
2657                         continue;
2658                 swap_div(bmap, i, pos);
2659                 --i;
2660         }
2661         return bmap;
2662 }
2663
2664 static int find_div(struct isl_basic_map *dst,
2665                         struct isl_basic_map *src, unsigned div)
2666 {
2667         int i;
2668
2669         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
2670
2671         for (i = 0; i < dst->n_div; ++i)
2672                 if (isl_seq_eq(dst->div[i], src->div[div], 1+1+total) &&
2673                     isl_seq_first_non_zero(dst->div[i]+1+1+total,
2674                                                 dst->n_div - src->n_div) == -1)
2675                         return i;
2676         return -1;
2677 }
2678
2679 /* For a div d = floor(f/m), add the constraints
2680  *
2681  *              f - m d >= 0
2682  *              -(f-(n-1)) + m d >= 0
2683  *
2684  * Note that the second constraint is the negation of
2685  *
2686  *              f - m d >= n
2687  */
2688 static int add_div_constraints(struct isl_ctx *ctx,
2689                 struct isl_basic_map *bmap, unsigned div)
2690 {
2691         int i, j;
2692         unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2693         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
2694
2695         i = isl_basic_map_alloc_inequality(ctx, bmap);
2696         if (i < 0)
2697                 return -1;
2698         isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
2699         isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
2700
2701         j = isl_basic_map_alloc_inequality(ctx, bmap);
2702         if (j < 0)
2703                 return -1;
2704         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
2705         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
2706         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
2707         return j;
2708 }
2709
2710 struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
2711                 struct isl_basic_map *dst, struct isl_basic_map *src)
2712 {
2713         int i;
2714         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
2715
2716         if (src->n_div == 0)
2717                 return dst;
2718
2719         src = order_divs(ctx, src);
2720         dst = isl_basic_map_extend(ctx, dst, dst->nparam, dst->n_in,
2721                         dst->n_out, src->n_div, 0, 2 * src->n_div);
2722         if (!dst)
2723                 return NULL;
2724         for (i = 0; i < src->n_div; ++i) {
2725                 int j = find_div(dst, src, i);
2726                 if (j < 0) {
2727                         j = isl_basic_map_alloc_div(ctx, dst);
2728                         if (j < 0)
2729                                 goto error;
2730                         isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
2731                         if (add_div_constraints(ctx, dst, j) < 0)
2732                                 goto error;
2733                 }
2734                 if (j != i)
2735                         swap_div(dst, i, j);
2736         }
2737         return dst;
2738 error:
2739         isl_basic_map_free(ctx, dst);
2740         return NULL;
2741 }
2742
2743 static struct isl_map *add_cut_constraint(struct isl_ctx *ctx,
2744                 struct isl_map *dst,
2745                 struct isl_basic_map *src, isl_int *c,
2746                 unsigned len, unsigned extra, int oppose)
2747 {
2748         struct isl_basic_map *copy = NULL;
2749         int is_empty;
2750         int k;
2751
2752         copy = isl_basic_map_copy(ctx, src);
2753         copy = isl_basic_map_cow(ctx, copy);
2754         if (!copy)
2755                 goto error;
2756         copy = isl_basic_map_extend(ctx, copy,
2757                 copy->nparam, copy->n_in, copy->n_out, 0, 0, 1);
2758         k = isl_basic_map_alloc_inequality(ctx, copy);
2759         if (k < 0)
2760                 goto error;
2761         if (oppose)
2762                 isl_seq_neg(copy->ineq[k], c, len);
2763         else
2764                 isl_seq_cpy(copy->ineq[k], c, len);
2765         isl_seq_clr(copy->ineq[k]+len, extra);
2766         isl_inequality_negate(ctx, copy, k);
2767         copy = isl_basic_map_simplify(ctx, copy);
2768         copy = isl_basic_map_finalize(ctx, copy);
2769         is_empty = isl_basic_map_is_empty(ctx, copy);
2770         if (is_empty < 0)
2771                 goto error;
2772         if (!is_empty)
2773                 dst = isl_map_add(ctx, dst, copy);
2774         else
2775                 isl_basic_map_free(ctx, copy);
2776         return dst;
2777 error:
2778         isl_basic_map_free(ctx, copy);
2779         isl_map_free(ctx, dst);
2780         return NULL;
2781 }
2782
2783 static struct isl_map *subtract(struct isl_ctx *ctx, struct isl_map *map,
2784                 struct isl_basic_map *bmap)
2785 {
2786         int i, j, k;
2787         unsigned flags = 0;
2788         struct isl_map *rest = NULL;
2789         unsigned max;
2790         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
2791
2792         assert(bmap);
2793
2794         if (!map)
2795                 goto error;
2796
2797         if (F_ISSET(map, ISL_MAP_DISJOINT))
2798                 FL_SET(flags, ISL_MAP_DISJOINT);
2799
2800         max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
2801         rest = isl_map_alloc(ctx, map->nparam, map->n_in, map->n_out,
2802                                 max, flags);
2803         if (!rest)
2804                 goto error;
2805
2806         for (i = 0; i < map->n; ++i) {
2807                 map->p[i] = isl_basic_map_align_divs(ctx, map->p[i], bmap);
2808                 if (!map->p[i])
2809                         goto error;
2810         }
2811
2812         for (j = 0; j < map->n; ++j)
2813                 map->p[j] = isl_basic_map_cow(ctx, map->p[j]);
2814
2815         for (i = 0; i < bmap->n_eq; ++i) {
2816                 for (j = 0; j < map->n; ++j) {
2817                         rest = add_cut_constraint(ctx, rest,
2818                                 map->p[j], bmap->eq[i],
2819                                 1+total, map->p[j]->n_div - bmap->n_div, 0);
2820                         if (!rest)
2821                                 goto error;
2822
2823                         rest = add_cut_constraint(ctx, rest,
2824                                 map->p[j], bmap->eq[i],
2825                                 1+total, map->p[j]->n_div - bmap->n_div, 1);
2826                         if (!rest)
2827                                 goto error;
2828
2829                         map->p[j] = isl_basic_map_extend(ctx, map->p[j],
2830                                 map->p[j]->nparam, map->p[j]->n_in,
2831                                 map->p[j]->n_out, 0, 1, 0);
2832                         if (!map->p[j])
2833                                 goto error;
2834                         k = isl_basic_map_alloc_equality(ctx, map->p[j]);
2835                         if (k < 0)
2836                                 goto error;
2837                         isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
2838                         isl_seq_clr(map->p[j]->eq[k]+1+total,
2839                                         map->p[j]->n_div - bmap->n_div);
2840                 }
2841         }
2842
2843         for (i = 0; i < bmap->n_ineq; ++i) {
2844                 for (j = 0; j < map->n; ++j) {
2845                         rest = add_cut_constraint(ctx, rest,
2846                                 map->p[j], bmap->ineq[i],
2847                                 1+total, map->p[j]->n_div - bmap->n_div, 0);
2848                         if (!rest)
2849                                 goto error;
2850
2851                         map->p[j] = isl_basic_map_extend(ctx, map->p[j],
2852                                 map->p[j]->nparam, map->p[j]->n_in,
2853                                 map->p[j]->n_out, 0, 0, 1);
2854                         if (!map->p[j])
2855                                 goto error;
2856                         k = isl_basic_map_alloc_inequality(ctx, map->p[j]);
2857                         if (k < 0)
2858                                 goto error;
2859                         isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
2860                         isl_seq_clr(map->p[j]->ineq[k]+1+total,
2861                                         map->p[j]->n_div - bmap->n_div);
2862                 }
2863         }
2864
2865         isl_map_free(ctx, map);
2866         return rest;
2867 error:
2868         isl_map_free(ctx, map);
2869         isl_map_free(ctx, rest);
2870         return NULL;
2871 }
2872
2873 struct isl_map *isl_map_subtract(struct isl_ctx *ctx, struct isl_map *map1,
2874                 struct isl_map *map2)
2875 {
2876         int i;
2877         if (!map1 || !map2)
2878                 goto error;
2879
2880         isl_assert(ctx, map1->nparam == map2->nparam, goto error);
2881         isl_assert(ctx, map1->n_in == map2->n_in, goto error);
2882         isl_assert(ctx, map1->n_out == map2->n_out, goto error);
2883
2884         if (isl_map_is_empty(ctx, map2)) {
2885                 isl_map_free(ctx, map2);
2886                 return map1;
2887         }
2888
2889         map1 = isl_map_compute_divs(ctx, map1);
2890         map2 = isl_map_compute_divs(ctx, map2);
2891         if (!map1 || !map2)
2892                 goto error;
2893
2894         for (i = 0; map1 && i < map2->n; ++i)
2895                 map1 = subtract(ctx, map1, map2->p[i]);
2896
2897         isl_map_free(ctx, map2);
2898         return map1;
2899 error:
2900         isl_map_free(ctx, map1);
2901         isl_map_free(ctx, map2);
2902         return NULL;
2903 }
2904
2905 struct isl_set *isl_set_subtract(struct isl_ctx *ctx, struct isl_set *set1,
2906                 struct isl_set *set2)
2907 {
2908         return (struct isl_set *)
2909                 isl_map_subtract(ctx,
2910                         (struct isl_map *)set1, (struct isl_map *)set2);
2911 }
2912
2913 struct isl_set *isl_set_apply(struct isl_ctx *ctx, struct isl_set *set,
2914                 struct isl_map *map)
2915 {
2916         isl_assert(ctx, set->dim == map->n_in, goto error);
2917         map = isl_map_intersect_domain(ctx, map, set);
2918         set = isl_map_range(ctx, map);
2919         return set;
2920 error:
2921         isl_set_free(ctx, set);
2922         isl_map_free(ctx, map);
2923         return NULL;
2924 }