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