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