isl_map.c

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