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