isl_map_simplify.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include "isl_equalities.h"
  11 #include "isl_map.h"
  12 #include "isl_map_private.h"
  13 #include "isl_seq.h"
  14 #include "isl_tab.h"
  15
  16 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
  17 {
  18         isl_int *t = bmap->eq[a];
  19         bmap->eq[a] = bmap->eq[b];
  20         bmap->eq[b] = t;
  21 }
  22
  23 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
  24 {
  25         if (a != b) {
  26                 isl_int *t = bmap->ineq[a];
  27                 bmap->ineq[a] = bmap->ineq[b];
  28                 bmap->ineq[b] = t;
  29         }
  30 }
  31
  32 static void set_swap_inequality(struct isl_basic_set *bset, int a, int b)
  33 {
  34         swap_inequality((struct isl_basic_map *)bset, a, b);
  35 }
  36
  37 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
  38 {
  39         isl_seq_cpy(c, c + n, rem);
  40         isl_seq_clr(c + rem, n);
  41 }
  42
  43 /* Drop n dimensions starting at first.
  44  *
  45  * In principle, this frees up some extra variables as the number
  46  * of columns remains constant, but we would have to extend
  47  * the div array too as the number of rows in this array is assumed
  48  * to be equal to extra.
  49  */
  50 struct isl_basic_set *isl_basic_set_drop_dims(
  51                 struct isl_basic_set *bset, unsigned first, unsigned n)
  52 {
  53         int i;
  54
  55         if (!bset)
  56                 goto error;
  57
  58         isl_assert(bset->ctx, first + n <= bset->dim->n_out, goto error);
  59
  60         if (n == 0)
  61                 return bset;
  62
  63         bset = isl_basic_set_cow(bset);
  64         if (!bset)
  65                 return NULL;
  66
  67         for (i = 0; i < bset->n_eq; ++i)
  68                 constraint_drop_vars(bset->eq[i]+1+bset->dim->nparam+first, n,
  69                                      (bset->dim->n_out-first-n)+bset->extra);
  70
  71         for (i = 0; i < bset->n_ineq; ++i)
  72                 constraint_drop_vars(bset->ineq[i]+1+bset->dim->nparam+first, n,
  73                                      (bset->dim->n_out-first-n)+bset->extra);
  74
  75         for (i = 0; i < bset->n_div; ++i)
  76                 constraint_drop_vars(bset->div[i]+1+1+bset->dim->nparam+first, n,
  77                                      (bset->dim->n_out-first-n)+bset->extra);
  78
  79         bset->dim = isl_dim_drop_outputs(bset->dim, first, n);
  80         if (!bset->dim)
  81                 goto error;
  82
  83         ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
  84         bset = isl_basic_set_simplify(bset);
  85         return isl_basic_set_finalize(bset);
  86 error:
  87         isl_basic_set_free(bset);
  88         return NULL;
  89 }
  90
  91 struct isl_set *isl_set_drop_dims(
  92                 struct isl_set *set, unsigned first, unsigned n)
  93 {
  94         int i;
  95
  96         if (!set)
  97                 goto error;
  98
  99         isl_assert(set->ctx, first + n <= set->dim->n_out, goto error);
 100
 101         if (n == 0)
 102                 return set;
 103         set = isl_set_cow(set);
 104         if (!set)
 105                 goto error;
 106         set->dim = isl_dim_drop_outputs(set->dim, first, n);
 107         if (!set->dim)
 108                 goto error;
 109
 110         for (i = 0; i < set->n; ++i) {
 111                 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
 112                 if (!set->p[i])
 113                         goto error;
 114         }
 115
 116         ISL_F_CLR(set, ISL_SET_NORMALIZED);
 117         return set;
 118 error:
 119         isl_set_free(set);
 120         return NULL;
 121 }
 122
 123 /* Move "n" divs starting at "first" to the end of the list of divs.
 124  */
 125 static struct isl_basic_map *move_divs_last(struct isl_basic_map *bmap,
 126         unsigned first, unsigned n)
 127 {
 128         isl_int **div;
 129         int i;
 130
 131         if (first + n == bmap->n_div)
 132                 return bmap;
 133
 134         div = isl_alloc_array(bmap->ctx, isl_int *, n);
 135         if (!div)
 136                 goto error;
 137         for (i = 0; i < n; ++i)
 138                 div[i] = bmap->div[first + i];
 139         for (i = 0; i < bmap->n_div - first - n; ++i)
 140                 bmap->div[first + i] = bmap->div[first + n + i];
 141         for (i = 0; i < n; ++i)
 142                 bmap->div[bmap->n_div - n + i] = div[i];
 143         free(div);
 144         return bmap;
 145 error:
 146         isl_basic_map_free(bmap);
 147         return NULL;
 148 }
 149
 150 /* Drop "n" dimensions of type "type" starting at "first".
 151  *
 152  * In principle, this frees up some extra variables as the number
 153  * of columns remains constant, but we would have to extend
 154  * the div array too as the number of rows in this array is assumed
 155  * to be equal to extra.
 156  */
 157 struct isl_basic_map *isl_basic_map_drop(struct isl_basic_map *bmap,
 158         enum isl_dim_type type, unsigned first, unsigned n)
 159 {
 160         int i;
 161         unsigned dim;
 162         unsigned offset;
 163         unsigned left;
 164
 165         if (!bmap)
 166                 goto error;
 167
 168         dim = isl_basic_map_dim(bmap, type);
 169         isl_assert(bmap->ctx, first + n <= dim, goto error);
 170
 171         if (n == 0)
 172                 return bmap;
 173
 174         bmap = isl_basic_map_cow(bmap);
 175         if (!bmap)
 176                 return NULL;
 177
 178         offset = isl_basic_map_offset(bmap, type) + first;
 179         left = isl_basic_map_total_dim(bmap) - (offset - 1) - n;
 180         for (i = 0; i < bmap->n_eq; ++i)
 181                 constraint_drop_vars(bmap->eq[i]+offset, n, left);
 182
 183         for (i = 0; i < bmap->n_ineq; ++i)
 184                 constraint_drop_vars(bmap->ineq[i]+offset, n, left);
 185
 186         for (i = 0; i < bmap->n_div; ++i)
 187                 constraint_drop_vars(bmap->div[i]+1+offset, n, left);
 188
 189         if (type == isl_dim_div) {
 190                 bmap = move_divs_last(bmap, first, n);
 191                 if (!bmap)
 192                         goto error;
 193                 isl_basic_map_free_div(bmap, n);
 194         } else
 195                 bmap->dim = isl_dim_drop(bmap->dim, type, first, n);
 196         if (!bmap->dim)
 197                 goto error;
 198
 199         ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
 200         bmap = isl_basic_map_simplify(bmap);
 201         return isl_basic_map_finalize(bmap);
 202 error:
 203         isl_basic_map_free(bmap);
 204         return NULL;
 205 }
 206
 207 struct isl_basic_map *isl_basic_map_drop_inputs(
 208                 struct isl_basic_map *bmap, unsigned first, unsigned n)
 209 {
 210         return isl_basic_map_drop(bmap, isl_dim_in, first, n);
 211 }
 212
 213 struct isl_map *isl_map_drop(struct isl_map *map,
 214         enum isl_dim_type type, unsigned first, unsigned n)
 215 {
 216         int i;
 217
 218         if (!map)
 219                 goto error;
 220
 221         isl_assert(map->ctx, first + n <= isl_map_dim(map, type), goto error);
 222
 223         if (n == 0)
 224                 return map;
 225         map = isl_map_cow(map);
 226         if (!map)
 227                 goto error;
 228         map->dim = isl_dim_drop(map->dim, type, first, n);
 229         if (!map->dim)
 230                 goto error;
 231
 232         for (i = 0; i < map->n; ++i) {
 233                 map->p[i] = isl_basic_map_drop(map->p[i], type, first, n);
 234                 if (!map->p[i])
 235                         goto error;
 236         }
 237         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
 238
 239         return map;
 240 error:
 241         isl_map_free(map);
 242         return NULL;
 243 }
 244
 245 struct isl_set *isl_set_drop(struct isl_set *set,
 246         enum isl_dim_type type, unsigned first, unsigned n)
 247 {
 248         return (isl_set *)isl_map_drop((isl_map *)set, type, first, n);
 249 }
 250
 251 struct isl_map *isl_map_drop_inputs(
 252                 struct isl_map *map, unsigned first, unsigned n)
 253 {
 254         return isl_map_drop(map, isl_dim_in, first, n);
 255 }
 256
 257 /*
 258  * We don't cow, as the div is assumed to be redundant.
 259  */
 260 static struct isl_basic_map *isl_basic_map_drop_div(
 261                 struct isl_basic_map *bmap, unsigned div)
 262 {
 263         int i;
 264         unsigned pos;
 265
 266         if (!bmap)
 267                 goto error;
 268
 269         pos = 1 + isl_dim_total(bmap->dim) + div;
 270
 271         isl_assert(bmap->ctx, div < bmap->n_div, goto error);
 272
 273         for (i = 0; i < bmap->n_eq; ++i)
 274                 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
 275
 276         for (i = 0; i < bmap->n_ineq; ++i) {
 277                 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
 278                         isl_basic_map_drop_inequality(bmap, i);
 279                         --i;
 280                         continue;
 281                 }
 282                 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
 283         }
 284
 285         for (i = 0; i < bmap->n_div; ++i)
 286                 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
 287
 288         if (div != bmap->n_div - 1) {
 289                 int j;
 290                 isl_int *t = bmap->div[div];
 291
 292                 for (j = div; j < bmap->n_div - 1; ++j)
 293                         bmap->div[j] = bmap->div[j+1];
 294
 295                 bmap->div[bmap->n_div - 1] = t;
 296         }
 297         ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
 298         isl_basic_map_free_div(bmap, 1);
 299
 300         return bmap;
 301 error:
 302         isl_basic_map_free(bmap);
 303         return NULL;
 304 }
 305
 306 struct isl_basic_map *isl_basic_map_normalize_constraints(
 307         struct isl_basic_map *bmap)
 308 {
 309         int i;
 310         isl_int gcd;
 311         unsigned total = isl_basic_map_total_dim(bmap);
 312
 313         isl_int_init(gcd);
 314         for (i = bmap->n_eq - 1; i >= 0; --i) {
 315                 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
 316                 if (isl_int_is_zero(gcd)) {
 317                         if (!isl_int_is_zero(bmap->eq[i][0])) {
 318                                 bmap = isl_basic_map_set_to_empty(bmap);
 319                                 break;
 320                         }
 321                         isl_basic_map_drop_equality(bmap, i);
 322                         continue;
 323                 }
 324                 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 325                         isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
 326                 if (isl_int_is_one(gcd))
 327                         continue;
 328                 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
 329                         bmap = isl_basic_map_set_to_empty(bmap);
 330                         break;
 331                 }
 332                 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
 333         }
 334
 335         for (i = bmap->n_ineq - 1; i >= 0; --i) {
 336                 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
 337                 if (isl_int_is_zero(gcd)) {
 338                         if (isl_int_is_neg(bmap->ineq[i][0])) {
 339                                 bmap = isl_basic_map_set_to_empty(bmap);
 340                                 break;
 341                         }
 342                         isl_basic_map_drop_inequality(bmap, i);
 343                         continue;
 344                 }
 345                 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 346                         isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
 347                 if (isl_int_is_one(gcd))
 348                         continue;
 349                 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
 350                 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
 351         }
 352         isl_int_clear(gcd);
 353
 354         return bmap;
 355 }
 356
 357 struct isl_basic_set *isl_basic_set_normalize_constraints(
 358         struct isl_basic_set *bset)
 359 {
 360         return (struct isl_basic_set *)isl_basic_map_normalize_constraints(
 361                 (struct isl_basic_map *)bset);
 362 }
 363
 364 /* Assumes divs have been ordered if keep_divs is set.
 365  */
 366 static void eliminate_var_using_equality(struct isl_basic_map *bmap,
 367         unsigned pos, isl_int *eq, int keep_divs, int *progress)
 368 {
 369         unsigned total;
 370         int k;
 371         int last_div;
 372
 373         total = isl_basic_map_total_dim(bmap);
 374         last_div = isl_seq_last_non_zero(eq + 1 + isl_dim_total(bmap->dim),
 375                                                 bmap->n_div);
 376         for (k = 0; k < bmap->n_eq; ++k) {
 377                 if (bmap->eq[k] == eq)
 378                         continue;
 379                 if (isl_int_is_zero(bmap->eq[k][1+pos]))
 380                         continue;
 381                 if (progress)
 382                         *progress = 1;
 383                 isl_seq_elim(bmap->eq[k], eq, 1+pos, 1+total, NULL);
 384         }
 385
 386         for (k = 0; k < bmap->n_ineq; ++k) {
 387                 if (isl_int_is_zero(bmap->ineq[k][1+pos]))
 388                         continue;
 389                 if (progress)
 390                         *progress = 1;
 391                 isl_seq_elim(bmap->ineq[k], eq, 1+pos, 1+total, NULL);
 392                 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
 393         }
 394
 395         for (k = 0; k < bmap->n_div; ++k) {
 396                 if (isl_int_is_zero(bmap->div[k][0]))
 397                         continue;
 398                 if (isl_int_is_zero(bmap->div[k][1+1+pos]))
 399                         continue;
 400                 if (progress)
 401                         *progress = 1;
 402                 /* We need to be careful about circular definitions,
 403                  * so for now we just remove the definition of div k
 404                  * if the equality contains any divs.
 405                  * If keep_divs is set, then the divs have been ordered
 406                  * and we can keep the definition as long as the result
 407                  * is still ordered.
 408                  */
 409                 if (last_div == -1 || (keep_divs && last_div < k))
 410                         isl_seq_elim(bmap->div[k]+1, eq,
 411                                         1+pos, 1+total, &bmap->div[k][0]);
 412                 else
 413                         isl_seq_clr(bmap->div[k], 1 + total);
 414                 ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
 415         }
 416 }
 417
 418 /* Assumes divs have been ordered if keep_divs is set.
 419  */
 420 static void eliminate_div(struct isl_basic_map *bmap, isl_int *eq,
 421         unsigned div, int keep_divs)
 422 {
 423         unsigned pos = isl_dim_total(bmap->dim) + div;
 424
 425         eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
 426
 427         isl_basic_map_drop_div(bmap, div);
 428 }
 429
 430 /* Elimininate divs based on equalities
 431  */
 432 static struct isl_basic_map *eliminate_divs_eq(
 433                 struct isl_basic_map *bmap, int *progress)
 434 {
 435         int d;
 436         int i;
 437         int modified = 0;
 438         unsigned off;
 439
 440         bmap = isl_basic_map_order_divs(bmap);
 441
 442         if (!bmap)
 443                 return NULL;
 444
 445         off = 1 + isl_dim_total(bmap->dim);
 446
 447         for (d = bmap->n_div - 1; d >= 0 ; --d) {
 448                 for (i = 0; i < bmap->n_eq; ++i) {
 449                         if (!isl_int_is_one(bmap->eq[i][off + d]) &&
 450                             !isl_int_is_negone(bmap->eq[i][off + d]))
 451                                 continue;
 452                         modified = 1;
 453                         *progress = 1;
 454                         eliminate_div(bmap, bmap->eq[i], d, 1);
 455                         isl_basic_map_drop_equality(bmap, i);
 456                         break;
 457                 }
 458         }
 459         if (modified)
 460                 return eliminate_divs_eq(bmap, progress);
 461         return bmap;
 462 }
 463
 464 /* Elimininate divs based on inequalities
 465  */
 466 static struct isl_basic_map *eliminate_divs_ineq(
 467                 struct isl_basic_map *bmap, int *progress)
 468 {
 469         int d;
 470         int i;
 471         unsigned off;
 472         struct isl_ctx *ctx;
 473
 474         if (!bmap)
 475                 return NULL;
 476
 477         ctx = bmap->ctx;
 478         off = 1 + isl_dim_total(bmap->dim);
 479
 480         for (d = bmap->n_div - 1; d >= 0 ; --d) {
 481                 for (i = 0; i < bmap->n_eq; ++i)
 482                         if (!isl_int_is_zero(bmap->eq[i][off + d]))
 483                                 break;
 484                 if (i < bmap->n_eq)
 485                         continue;
 486                 for (i = 0; i < bmap->n_ineq; ++i)
 487                         if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
 488                                 break;
 489                 if (i < bmap->n_ineq)
 490                         continue;
 491                 *progress = 1;
 492                 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
 493                 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
 494                         break;
 495                 bmap = isl_basic_map_drop_div(bmap, d);
 496                 if (!bmap)
 497                         break;
 498         }
 499         return bmap;
 500 }
 501
 502 struct isl_basic_map *isl_basic_map_gauss(
 503         struct isl_basic_map *bmap, int *progress)
 504 {
 505         int k;
 506         int done;
 507         int last_var;
 508         unsigned total_var;
 509         unsigned total;
 510
 511         bmap = isl_basic_map_order_divs(bmap);
 512
 513         if (!bmap)
 514                 return NULL;
 515
 516         total = isl_basic_map_total_dim(bmap);
 517         total_var = total - bmap->n_div;
 518
 519         last_var = total - 1;
 520         for (done = 0; done < bmap->n_eq; ++done) {
 521                 for (; last_var >= 0; --last_var) {
 522                         for (k = done; k < bmap->n_eq; ++k)
 523                                 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
 524                                         break;
 525                         if (k < bmap->n_eq)
 526                                 break;
 527                 }
 528                 if (last_var < 0)
 529                         break;
 530                 if (k != done)
 531                         swap_equality(bmap, k, done);
 532                 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
 533                         isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
 534
 535                 eliminate_var_using_equality(bmap, last_var, bmap->eq[done], 1,
 536                                                 progress);
 537
 538                 if (last_var >= total_var &&
 539                     isl_int_is_zero(bmap->div[last_var - total_var][0])) {
 540                         unsigned div = last_var - total_var;
 541                         isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
 542                         isl_int_set_si(bmap->div[div][1+1+last_var], 0);
 543                         isl_int_set(bmap->div[div][0],
 544                                     bmap->eq[done][1+last_var]);
 545                         ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
 546                 }
 547         }
 548         if (done == bmap->n_eq)
 549                 return bmap;
 550         for (k = done; k < bmap->n_eq; ++k) {
 551                 if (isl_int_is_zero(bmap->eq[k][0]))
 552                         continue;
 553                 return isl_basic_map_set_to_empty(bmap);
 554         }
 555         isl_basic_map_free_equality(bmap, bmap->n_eq-done);
 556         return bmap;
 557 }
 558
 559 struct isl_basic_set *isl_basic_set_gauss(
 560         struct isl_basic_set *bset, int *progress)
 561 {
 562         return (struct isl_basic_set*)isl_basic_map_gauss(
 563                         (struct isl_basic_map *)bset, progress);
 564 }
 565
 566
 567 static unsigned int round_up(unsigned int v)
 568 {
 569         int old_v = v;
 570
 571         while (v) {
 572                 old_v = v;
 573                 v ^= v & -v;
 574         }
 575         return old_v << 1;
 576 }
 577
 578 static int hash_index(isl_int ***index, unsigned int size, int bits,
 579                         struct isl_basic_map *bmap, int k)
 580 {
 581         int h;
 582         unsigned total = isl_basic_map_total_dim(bmap);
 583         uint32_t hash = isl_seq_get_hash_bits(bmap->ineq[k]+1, total, bits);
 584         for (h = hash; index[h]; h = (h+1) % size)
 585                 if (&bmap->ineq[k] != index[h] &&
 586                     isl_seq_eq(bmap->ineq[k]+1, index[h][0]+1, total))
 587                         break;
 588         return h;
 589 }
 590
 591 static int set_hash_index(isl_int ***index, unsigned int size, int bits,
 592                           struct isl_basic_set *bset, int k)
 593 {
 594         return hash_index(index, size, bits, (struct isl_basic_map *)bset, k);
 595 }
 596
 597 /* If we can eliminate more than one div, then we need to make
 598  * sure we do it from last div to first div, in order not to
 599  * change the position of the other divs that still need to
 600  * be removed.
 601  */
 602 static struct isl_basic_map *remove_duplicate_divs(
 603         struct isl_basic_map *bmap, int *progress)
 604 {
 605         unsigned int size;
 606         int *index;
 607         int *elim_for;
 608         int k, l, h;
 609         int bits;
 610         struct isl_blk eq;
 611         unsigned total_var = isl_dim_total(bmap->dim);
 612         unsigned total = total_var + bmap->n_div;
 613         struct isl_ctx *ctx;
 614
 615         if (bmap->n_div <= 1)
 616                 return bmap;
 617
 618         ctx = bmap->ctx;
 619         for (k = bmap->n_div - 1; k >= 0; --k)
 620                 if (!isl_int_is_zero(bmap->div[k][0]))
 621                         break;
 622         if (k <= 0)
 623                 return bmap;
 624
 625         elim_for = isl_calloc_array(ctx, int, bmap->n_div);
 626         size = round_up(4 * bmap->n_div / 3 - 1);
 627         bits = ffs(size) - 1;
 628         index = isl_calloc_array(ctx, int, size);
 629         if (!index)
 630                 return bmap;
 631         eq = isl_blk_alloc(ctx, 1+total);
 632         if (isl_blk_is_error(eq))
 633                 goto out;
 634
 635         isl_seq_clr(eq.data, 1+total);
 636         index[isl_seq_get_hash_bits(bmap->div[k], 2+total, bits)] = k + 1;
 637         for (--k; k >= 0; --k) {
 638                 uint32_t hash;
 639
 640                 if (isl_int_is_zero(bmap->div[k][0]))
 641                         continue;
 642
 643                 hash = isl_seq_get_hash_bits(bmap->div[k], 2+total, bits);
 644                 for (h = hash; index[h]; h = (h+1) % size)
 645                         if (isl_seq_eq(bmap->div[k],
 646                                        bmap->div[index[h]-1], 2+total))
 647                                 break;
 648                 if (index[h]) {
 649                         *progress = 1;
 650                         l = index[h] - 1;
 651                         elim_for[l] = k + 1;
 652                 }
 653                 index[h] = k+1;
 654         }
 655         for (l = bmap->n_div - 1; l >= 0; --l) {
 656                 if (!elim_for[l])
 657                         continue;
 658                 k = elim_for[l] - 1;
 659                 isl_int_set_si(eq.data[1+total_var+k], -1);
 660                 isl_int_set_si(eq.data[1+total_var+l], 1);
 661                 eliminate_div(bmap, eq.data, l, 0);
 662                 isl_int_set_si(eq.data[1+total_var+k], 0);
 663                 isl_int_set_si(eq.data[1+total_var+l], 0);
 664         }
 665
 666         isl_blk_free(ctx, eq);
 667 out:
 668         free(index);
 669         free(elim_for);
 670         return bmap;
 671 }
 672
 673 static int n_pure_div_eq(struct isl_basic_map *bmap)
 674 {
 675         int i, j;
 676         unsigned total;
 677
 678         total = isl_dim_total(bmap->dim);
 679         for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
 680                 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
 681                         --j;
 682                 if (j < 0)
 683                         break;
 684                 if (isl_seq_first_non_zero(bmap->eq[i] + 1 + total, j) != -1)
 685                         return 0;
 686         }
 687         return i;
 688 }
 689
 690 /* Normalize divs that appear in equalities.
 691  *
 692  * In particular, we assume that bmap contains some equalities
 693  * of the form
 694  *
 695  *      a x = m * e_i
 696  *
 697  * and we want to replace the set of e_i by a minimal set and
 698  * such that the new e_i have a canonical representation in terms
 699  * of the vector x.
 700  * If any of the equalities involves more than one divs, then
 701  * we currently simply bail out.
 702  *
 703  * Let us first additionally assume that all equalities involve
 704  * a div.  The equalities then express modulo constraints on the
 705  * remaining variables and we can use "parameter compression"
 706  * to find a minimal set of constraints.  The result is a transformation
 707  *
 708  *      x = T(x') = x_0 + G x'
 709  *
 710  * with G a lower-triangular matrix with all elements below the diagonal
 711  * non-negative and smaller than the diagonal element on the same row.
 712  * We first normalize x_0 by making the same property hold in the affine
 713  * T matrix.
 714  * The rows i of G with a 1 on the diagonal do not impose any modulo
 715  * constraint and simply express x_i = x'_i.
 716  * For each of the remaining rows i, we introduce a div and a corresponding
 717  * equality.  In particular
 718  *
 719  *      g_ii e_j = x_i - g_i(x')
 720  *
 721  * where each x'_k is replaced either by x_k (if g_kk = 1) or the
 722  * corresponding div (if g_kk != 1).
 723  *
 724  * If there are any equalities not involving any div, then we
 725  * first apply a variable compression on the variables x:
 726  *
 727  *      x = C x''       x'' = C_2 x
 728  *
 729  * and perform the above parameter compression on A C instead of on A.
 730  * The resulting compression is then of the form
 731  *
 732  *      x'' = T(x') = x_0 + G x'
 733  *
 734  * and in constructing the new divs and the corresponding equalities,
 735  * we have to replace each x'', i.e., the x'_k with (g_kk = 1),
 736  * by the corresponding row from C_2.
 737  */
 738 static struct isl_basic_map *normalize_divs(
 739         struct isl_basic_map *bmap, int *progress)
 740 {
 741         int i, j, k;
 742         int total;
 743         int div_eq;
 744         struct isl_mat *B;
 745         struct isl_vec *d;
 746         struct isl_mat *T = NULL;
 747         struct isl_mat *C = NULL;
 748         struct isl_mat *C2 = NULL;
 749         isl_int v;
 750         int *pos;
 751         int dropped, needed;
 752
 753         if (!bmap)
 754                 return NULL;
 755
 756         if (bmap->n_div == 0)
 757                 return bmap;
 758
 759         if (bmap->n_eq == 0)
 760                 return bmap;
 761
 762         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
 763                 return bmap;
 764
 765         total = isl_dim_total(bmap->dim);
 766         div_eq = n_pure_div_eq(bmap);
 767         if (div_eq == 0)
 768                 return bmap;
 769
 770         if (div_eq < bmap->n_eq) {
 771                 B = isl_mat_sub_alloc(bmap->ctx, bmap->eq, div_eq,
 772                                         bmap->n_eq - div_eq, 0, 1 + total);
 773                 C = isl_mat_variable_compression(B, &C2);
 774                 if (!C || !C2)
 775                         goto error;
 776                 if (C->n_col == 0) {
 777                         bmap = isl_basic_map_set_to_empty(bmap);
 778                         isl_mat_free(C);
 779                         isl_mat_free(C2);
 780                         goto done;
 781                 }
 782         }
 783
 784         d = isl_vec_alloc(bmap->ctx, div_eq);
 785         if (!d)
 786                 goto error;
 787         for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
 788                 while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + total + j]))
 789                         --j;
 790                 isl_int_set(d->block.data[i], bmap->eq[i][1 + total + j]);
 791         }
 792         B = isl_mat_sub_alloc(bmap->ctx, bmap->eq, 0, div_eq, 0, 1 + total);
 793
 794         if (C) {
 795                 B = isl_mat_product(B, C);
 796                 C = NULL;
 797         }
 798
 799         T = isl_mat_parameter_compression(B, d);
 800         if (!T)
 801                 goto error;
 802         if (T->n_col == 0) {
 803                 bmap = isl_basic_map_set_to_empty(bmap);
 804                 isl_mat_free(C2);
 805                 isl_mat_free(T);
 806                 goto done;
 807         }
 808         isl_int_init(v);
 809         for (i = 0; i < T->n_row - 1; ++i) {
 810                 isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
 811                 if (isl_int_is_zero(v))
 812                         continue;
 813                 isl_mat_col_submul(T, 0, v, 1 + i);
 814         }
 815         isl_int_clear(v);
 816         pos = isl_alloc_array(bmap->ctx, int, T->n_row);
 817         /* We have to be careful because dropping equalities may reorder them */
 818         dropped = 0;
 819         for (j = bmap->n_div - 1; j >= 0; --j) {
 820                 for (i = 0; i < bmap->n_eq; ++i)
 821                         if (!isl_int_is_zero(bmap->eq[i][1 + total + j]))
 822                                 break;
 823                 if (i < bmap->n_eq) {
 824                         bmap = isl_basic_map_drop_div(bmap, j);
 825                         isl_basic_map_drop_equality(bmap, i);
 826                         ++dropped;
 827                 }
 828         }
 829         pos[0] = 0;
 830         needed = 0;
 831         for (i = 1; i < T->n_row; ++i) {
 832                 if (isl_int_is_one(T->row[i][i]))
 833                         pos[i] = i;
 834                 else
 835                         needed++;
 836         }
 837         if (needed > dropped) {
 838                 bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
 839                                 needed, needed, 0);
 840                 if (!bmap)
 841                         goto error;
 842         }
 843         for (i = 1; i < T->n_row; ++i) {
 844                 if (isl_int_is_one(T->row[i][i]))
 845                         continue;
 846                 k = isl_basic_map_alloc_div(bmap);
 847                 pos[i] = 1 + total + k;
 848                 isl_seq_clr(bmap->div[k] + 1, 1 + total + bmap->n_div);
 849                 isl_int_set(bmap->div[k][0], T->row[i][i]);
 850                 if (C2)
 851                         isl_seq_cpy(bmap->div[k] + 1, C2->row[i], 1 + total);
 852                 else
 853                         isl_int_set_si(bmap->div[k][1 + i], 1);
 854                 for (j = 0; j < i; ++j) {
 855                         if (isl_int_is_zero(T->row[i][j]))
 856                                 continue;
 857                         if (pos[j] < T->n_row && C2)
 858                                 isl_seq_submul(bmap->div[k] + 1, T->row[i][j],
 859                                                 C2->row[pos[j]], 1 + total);
 860                         else
 861                                 isl_int_neg(bmap->div[k][1 + pos[j]],
 862                                                                 T->row[i][j]);
 863                 }
 864                 j = isl_basic_map_alloc_equality(bmap);
 865                 isl_seq_neg(bmap->eq[j], bmap->div[k]+1, 1+total+bmap->n_div);
 866                 isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
 867         }
 868         free(pos);
 869         isl_mat_free(C2);
 870         isl_mat_free(T);
 871
 872         if (progress)
 873                 *progress = 1;
 874 done:
 875         ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
 876
 877         return bmap;
 878 error:
 879         isl_mat_free(C);
 880         isl_mat_free(C2);
 881         isl_mat_free(T);
 882         return bmap;
 883 }
 884
 885 static struct isl_basic_map *set_div_from_lower_bound(
 886         struct isl_basic_map *bmap, int div, int ineq)
 887 {
 888         unsigned total = 1 + isl_dim_total(bmap->dim);
 889
 890         isl_seq_neg(bmap->div[div] + 1, bmap->ineq[ineq], total + bmap->n_div);
 891         isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
 892         isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
 893         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
 894         isl_int_set_si(bmap->div[div][1 + total + div], 0);
 895
 896         return bmap;
 897 }
 898
 899 /* Check whether it is ok to define a div based on an inequality.
 900  * To avoid the introduction of circular definitions of divs, we
 901  * do not allow such a definition if the resulting expression would refer to
 902  * any other undefined divs or if any known div is defined in
 903  * terms of the unknown div.
 904  */
 905 static int ok_to_set_div_from_bound(struct isl_basic_map *bmap,
 906         int div, int ineq)
 907 {
 908         int j;
 909         unsigned total = 1 + isl_dim_total(bmap->dim);
 910
 911         /* Not defined in terms of unknown divs */
 912         for (j = 0; j < bmap->n_div; ++j) {
 913                 if (div == j)
 914                         continue;
 915                 if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
 916                         continue;
 917                 if (isl_int_is_zero(bmap->div[j][0]))
 918                         return 0;
 919         }
 920
 921         /* No other div defined in terms of this one => avoid loops */
 922         for (j = 0; j < bmap->n_div; ++j) {
 923                 if (div == j)
 924                         continue;
 925                 if (isl_int_is_zero(bmap->div[j][0]))
 926                         continue;
 927                 if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
 928                         return 0;
 929         }
 930
 931         return 1;
 932 }
 933
 934 /* Given two constraints "k" and "l" that are opposite to each other,
 935  * except for the constant term, check if we can use them
 936  * to obtain an expression for one of the hitherto unknown divs.
 937  * "sum" is the sum of the constant terms of the constraints.
 938  * If this sum is strictly smaller than the coefficient of one
 939  * of the divs, then this pair can be used define the div.
 940  * To avoid the introduction of circular definitions of divs, we
 941  * do not use the pair if the resulting expression would refer to
 942  * any other undefined divs or if any known div is defined in
 943  * terms of the unknown div.
 944  */
 945 static struct isl_basic_map *check_for_div_constraints(
 946         struct isl_basic_map *bmap, int k, int l, isl_int sum, int *progress)
 947 {
 948         int i;
 949         unsigned total = 1 + isl_dim_total(bmap->dim);
 950
 951         for (i = 0; i < bmap->n_div; ++i) {
 952                 if (!isl_int_is_zero(bmap->div[i][0]))
 953                         continue;
 954                 if (isl_int_is_zero(bmap->ineq[k][total + i]))
 955                         continue;
 956                 if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
 957                         continue;
 958                 if (!ok_to_set_div_from_bound(bmap, i, k))
 959                         break;
 960                 if (isl_int_is_pos(bmap->ineq[k][total + i]))
 961                         bmap = set_div_from_lower_bound(bmap, i, k);
 962                 else
 963                         bmap = set_div_from_lower_bound(bmap, i, l);
 964                 if (progress)
 965                         *progress = 1;
 966                 break;
 967         }
 968         return bmap;
 969 }
 970
 971 static struct isl_basic_map *remove_duplicate_constraints(
 972         struct isl_basic_map *bmap, int *progress)
 973 {
 974         unsigned int size;
 975         isl_int ***index;
 976         int k, l, h;
 977         int bits;
 978         unsigned total = isl_basic_map_total_dim(bmap);
 979         isl_int sum;
 980
 981         if (bmap->n_ineq <= 1)
 982                 return bmap;
 983
 984         size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
 985         bits = ffs(size) - 1;
 986         index = isl_calloc_array(ctx, isl_int **, size);
 987         if (!index)
 988                 return bmap;
 989
 990         index[isl_seq_get_hash_bits(bmap->ineq[0]+1, total, bits)] = &bmap->ineq[0];
 991         for (k = 1; k < bmap->n_ineq; ++k) {
 992                 h = hash_index(index, size, bits, bmap, k);
 993                 if (!index[h]) {
 994                         index[h] = &bmap->ineq[k];
 995                         continue;
 996                 }
 997                 if (progress)
 998                         *progress = 1;
 999                 l = index[h] - &bmap->ineq[0];
1000                 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1001                         swap_inequality(bmap, k, l);
1002                 isl_basic_map_drop_inequality(bmap, k);
1003                 --k;
1004         }
1005         isl_int_init(sum);
1006         for (k = 0; k < bmap->n_ineq-1; ++k) {
1007                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1008                 h = hash_index(index, size, bits, bmap, k);
1009                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
1010                 if (!index[h])
1011                         continue;
1012                 l = index[h] - &bmap->ineq[0];
1013                 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1014                 if (isl_int_is_pos(sum)) {
1015                         bmap = check_for_div_constraints(bmap, k, l, sum,
1016                                                          progress);
1017                         continue;
1018                 }
1019                 if (isl_int_is_zero(sum)) {
1020                         /* We need to break out of the loop after these
1021                          * changes since the contents of the hash
1022                          * will no longer be valid.
1023                          * Plus, we probably we want to regauss first.
1024                          */
1025                         if (progress)
1026                                 *progress = 1;
1027                         isl_basic_map_drop_inequality(bmap, l);
1028                         isl_basic_map_inequality_to_equality(bmap, k);
1029                 } else
1030                         bmap = isl_basic_map_set_to_empty(bmap);
1031                 break;
1032         }
1033         isl_int_clear(sum);
1034
1035         free(index);
1036         return bmap;
1037 }
1038
1039
1040 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1041 {
1042         int progress = 1;
1043         if (!bmap)
1044                 return NULL;
1045         while (progress) {
1046                 progress = 0;
1047                 bmap = isl_basic_map_normalize_constraints(bmap);
1048                 bmap = remove_duplicate_divs(bmap, &progress);
1049                 bmap = eliminate_divs_eq(bmap, &progress);
1050                 bmap = eliminate_divs_ineq(bmap, &progress);
1051                 bmap = isl_basic_map_gauss(bmap, &progress);
1052                 /* requires equalities in normal form */
1053                 bmap = normalize_divs(bmap, &progress);
1054                 bmap = remove_duplicate_constraints(bmap, &progress);
1055         }
1056         return bmap;
1057 }
1058
1059 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1060 {
1061         return (struct isl_basic_set *)
1062                 isl_basic_map_simplify((struct isl_basic_map *)bset);
1063 }
1064
1065
1066 /* If the only constraints a div d=floor(f/m)
1067  * appears in are its two defining constraints
1068  *
1069  *      f - m d >=0
1070  *      -(f - (m - 1)) + m d >= 0
1071  *
1072  * then it can safely be removed.
1073  */
1074 static int div_is_redundant(struct isl_basic_map *bmap, int div)
1075 {
1076         int i;
1077         unsigned pos = 1 + isl_dim_total(bmap->dim) + div;
1078
1079         for (i = 0; i < bmap->n_eq; ++i)
1080                 if (!isl_int_is_zero(bmap->eq[i][pos]))
1081                         return 0;
1082
1083         for (i = 0; i < bmap->n_ineq; ++i) {
1084                 if (isl_int_is_zero(bmap->ineq[i][pos]))
1085                         continue;
1086                 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
1087                         int neg;
1088                         isl_int_sub(bmap->div[div][1],
1089                                         bmap->div[div][1], bmap->div[div][0]);
1090                         isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1091                         neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
1092                         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1093                         isl_int_add(bmap->div[div][1],
1094                                         bmap->div[div][1], bmap->div[div][0]);
1095                         if (!neg)
1096                                 return 0;
1097                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
1098                                                     bmap->n_div-div-1) != -1)
1099                                 return 0;
1100                 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
1101                         if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
1102                                 return 0;
1103                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
1104                                                     bmap->n_div-div-1) != -1)
1105                                 return 0;
1106                 } else
1107                         return 0;
1108         }
1109
1110         for (i = 0; i < bmap->n_div; ++i)
1111                 if (!isl_int_is_zero(bmap->div[i][1+pos]))
1112                         return 0;
1113
1114         return 1;
1115 }
1116
1117 /*
1118  * Remove divs that don't occur in any of the constraints or other divs.
1119  * These can arise when dropping some of the variables in a quast
1120  * returned by piplib.
1121  */
1122 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
1123 {
1124         int i;
1125
1126         if (!bmap)
1127                 return NULL;
1128
1129         for (i = bmap->n_div-1; i >= 0; --i) {
1130                 if (!div_is_redundant(bmap, i))
1131                         continue;
1132                 bmap = isl_basic_map_drop_div(bmap, i);
1133         }
1134         return bmap;
1135 }
1136
1137 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
1138 {
1139         bmap = remove_redundant_divs(bmap);
1140         if (!bmap)
1141                 return NULL;
1142         ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1143         return bmap;
1144 }
1145
1146 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
1147 {
1148         return (struct isl_basic_set *)
1149                 isl_basic_map_finalize((struct isl_basic_map *)bset);
1150 }
1151
1152 struct isl_set *isl_set_finalize(struct isl_set *set)
1153 {
1154         int i;
1155
1156         if (!set)
1157                 return NULL;
1158         for (i = 0; i < set->n; ++i) {
1159                 set->p[i] = isl_basic_set_finalize(set->p[i]);
1160                 if (!set->p[i])
1161                         goto error;
1162         }
1163         return set;
1164 error:
1165         isl_set_free(set);
1166         return NULL;
1167 }
1168
1169 struct isl_map *isl_map_finalize(struct isl_map *map)
1170 {
1171         int i;
1172
1173         if (!map)
1174                 return NULL;
1175         for (i = 0; i < map->n; ++i) {
1176                 map->p[i] = isl_basic_map_finalize(map->p[i]);
1177                 if (!map->p[i])
1178                         goto error;
1179         }
1180         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1181         return map;
1182 error:
1183         isl_map_free(map);
1184         return NULL;
1185 }
1186
1187
1188 /* Remove definition of any div that is defined in terms of the given variable.
1189  * The div itself is not removed.  Functions such as
1190  * eliminate_divs_ineq depend on the other divs remaining in place.
1191  */
1192 static struct isl_basic_map *remove_dependent_vars(struct isl_basic_map *bmap,
1193                                                                         int pos)
1194 {
1195         int i;
1196
1197         for (i = 0; i < bmap->n_div; ++i) {
1198                 if (isl_int_is_zero(bmap->div[i][0]))
1199                         continue;
1200                 if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1201                         continue;
1202                 isl_int_set_si(bmap->div[i][0], 0);
1203         }
1204         return bmap;
1205 }
1206
1207 /* Eliminate the specified variables from the constraints using
1208  * Fourier-Motzkin.  The variables themselves are not removed.
1209  */
1210 struct isl_basic_map *isl_basic_map_eliminate_vars(
1211         struct isl_basic_map *bmap, unsigned pos, unsigned n)
1212 {
1213         int d;
1214         int i, j, k;
1215         unsigned total;
1216
1217         if (n == 0)
1218                 return bmap;
1219         if (!bmap)
1220                 return NULL;
1221         total = isl_basic_map_total_dim(bmap);
1222
1223         bmap = isl_basic_map_cow(bmap);
1224         for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1225                 bmap = remove_dependent_vars(bmap, d);
1226
1227         for (d = pos + n - 1;
1228              d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1229                 isl_seq_clr(bmap->div[d-(total-bmap->n_div)], 2+total);
1230         for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1231                 int n_lower, n_upper;
1232                 if (!bmap)
1233                         return NULL;
1234                 for (i = 0; i < bmap->n_eq; ++i) {
1235                         if (isl_int_is_zero(bmap->eq[i][1+d]))
1236                                 continue;
1237                         eliminate_var_using_equality(bmap, d, bmap->eq[i], 0, NULL);
1238                         isl_basic_map_drop_equality(bmap, i);
1239                         break;
1240                 }
1241                 if (i < bmap->n_eq)
1242                         continue;
1243                 n_lower = 0;
1244                 n_upper = 0;
1245                 for (i = 0; i < bmap->n_ineq; ++i) {
1246                         if (isl_int_is_pos(bmap->ineq[i][1+d]))
1247                                 n_lower++;
1248                         else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1249                                 n_upper++;
1250                 }
1251                 bmap = isl_basic_map_extend_constraints(bmap,
1252                                 0, n_lower * n_upper);
1253                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1254                         int last;
1255                         if (isl_int_is_zero(bmap->ineq[i][1+d]))
1256                                 continue;
1257                         last = -1;
1258                         for (j = 0; j < i; ++j) {
1259                                 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1260                                         continue;
1261                                 last = j;
1262                                 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1263                                     isl_int_sgn(bmap->ineq[j][1+d]))
1264                                         continue;
1265                                 k = isl_basic_map_alloc_inequality(bmap);
1266                                 if (k < 0)
1267                                         goto error;
1268                                 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1269                                                 1+total);
1270                                 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1271                                                 1+d, 1+total, NULL);
1272                         }
1273                         isl_basic_map_drop_inequality(bmap, i);
1274                         i = last + 1;
1275                 }
1276                 if (n_lower > 0 && n_upper > 0) {
1277                         bmap = isl_basic_map_normalize_constraints(bmap);
1278                         bmap = remove_duplicate_constraints(bmap, NULL);
1279                         bmap = isl_basic_map_gauss(bmap, NULL);
1280                         bmap = isl_basic_map_convex_hull(bmap);
1281                         if (!bmap)
1282                                 goto error;
1283                         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1284                                 break;
1285                 }
1286         }
1287         ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED);
1288         return bmap;
1289 error:
1290         isl_basic_map_free(bmap);
1291         return NULL;
1292 }
1293
1294 struct isl_basic_set *isl_basic_set_eliminate_vars(
1295         struct isl_basic_set *bset, unsigned pos, unsigned n)
1296 {
1297         return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1298                         (struct isl_basic_map *)bset, pos, n);
1299 }
1300
1301 /* Don't assume equalities are in order, because align_divs
1302  * may have changed the order of the divs.
1303  */
1304 static void compute_elimination_index(struct isl_basic_map *bmap, int *elim)
1305 {
1306         int d, i;
1307         unsigned total;
1308
1309         total = isl_dim_total(bmap->dim);
1310         for (d = 0; d < total; ++d)
1311                 elim[d] = -1;
1312         for (i = 0; i < bmap->n_eq; ++i) {
1313                 for (d = total - 1; d >= 0; --d) {
1314                         if (isl_int_is_zero(bmap->eq[i][1+d]))
1315                                 continue;
1316                         elim[d] = i;
1317                         break;
1318                 }
1319         }
1320 }
1321
1322 static void set_compute_elimination_index(struct isl_basic_set *bset, int *elim)
1323 {
1324         compute_elimination_index((struct isl_basic_map *)bset, elim);
1325 }
1326
1327 static int reduced_using_equalities(isl_int *dst, isl_int *src,
1328         struct isl_basic_map *bmap, int *elim)
1329 {
1330         int d;
1331         int copied = 0;
1332         unsigned total;
1333
1334         total = isl_dim_total(bmap->dim);
1335         for (d = total - 1; d >= 0; --d) {
1336                 if (isl_int_is_zero(src[1+d]))
1337                         continue;
1338                 if (elim[d] == -1)
1339                         continue;
1340                 if (!copied) {
1341                         isl_seq_cpy(dst, src, 1 + total);
1342                         copied = 1;
1343                 }
1344                 isl_seq_elim(dst, bmap->eq[elim[d]], 1 + d, 1 + total, NULL);
1345         }
1346         return copied;
1347 }
1348
1349 static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
1350         struct isl_basic_set *bset, int *elim)
1351 {
1352         return reduced_using_equalities(dst, src,
1353                                         (struct isl_basic_map *)bset, elim);
1354 }
1355
1356 static struct isl_basic_set *isl_basic_set_reduce_using_equalities(
1357         struct isl_basic_set *bset, struct isl_basic_set *context)
1358 {
1359         int i;
1360         int *elim;
1361
1362         if (!bset || !context)
1363                 goto error;
1364
1365         bset = isl_basic_set_cow(bset);
1366         if (!bset)
1367                 goto error;
1368
1369         elim = isl_alloc_array(ctx, int, isl_basic_set_n_dim(bset));
1370         if (!elim)
1371                 goto error;
1372         set_compute_elimination_index(context, elim);
1373         for (i = 0; i < bset->n_eq; ++i)
1374                 set_reduced_using_equalities(bset->eq[i], bset->eq[i],
1375                                                         context, elim);
1376         for (i = 0; i < bset->n_ineq; ++i)
1377                 set_reduced_using_equalities(bset->ineq[i], bset->ineq[i],
1378                                                         context, elim);
1379         isl_basic_set_free(context);
1380         free(elim);
1381         bset = isl_basic_set_simplify(bset);
1382         bset = isl_basic_set_finalize(bset);
1383         return bset;
1384 error:
1385         isl_basic_set_free(bset);
1386         isl_basic_set_free(context);
1387         return NULL;
1388 }
1389
1390 static struct isl_basic_set *remove_shifted_constraints(
1391         struct isl_basic_set *bset, struct isl_basic_set *context)
1392 {
1393         unsigned int size;
1394         isl_int ***index;
1395         int bits;
1396         int k, h, l;
1397
1398         if (!bset)
1399                 return NULL;
1400
1401         size = round_up(4 * (context->n_ineq+1) / 3 - 1);
1402         bits = ffs(size) - 1;
1403         index = isl_calloc_array(ctx, isl_int **, size);
1404         if (!index)
1405                 return bset;
1406
1407         for (k = 0; k < context->n_ineq; ++k) {
1408                 h = set_hash_index(index, size, bits, context, k);
1409                 index[h] = &context->ineq[k];
1410         }
1411         for (k = 0; k < bset->n_ineq; ++k) {
1412                 h = set_hash_index(index, size, bits, bset, k);
1413                 if (!index[h])
1414                         continue;
1415                 l = index[h] - &context->ineq[0];
1416                 if (isl_int_lt(bset->ineq[k][0], context->ineq[l][0]))
1417                         continue;
1418                 bset = isl_basic_set_cow(bset);
1419                 if (!bset)
1420                         goto error;
1421                 isl_basic_set_drop_inequality(bset, k);
1422                 --k;
1423         }
1424         free(index);
1425         return bset;
1426 error:
1427         free(index);
1428         return bset;
1429 }
1430
1431 /* Tighten (decrease) the constant terms of the inequalities based
1432  * on the equalities, without removing any integer points.
1433  * For example, if there is an equality
1434  *
1435  *              i = 3 * j
1436  *
1437  * and an inequality
1438  *
1439  *              i >= 1
1440  *
1441  * then we want to replace the inequality by
1442  *
1443  *              i >= 3
1444  *
1445  * We do this by computing a variable compression and translating
1446  * the constraints to the compressed space.
1447  * If any constraint has coefficients (except the contant term)
1448  * with a common factor "f", then we can replace the constant term "c"
1449  * by
1450  *
1451  *              f * floor(c/f)
1452  *
1453  * That is, we add
1454  *
1455  *              f * floor(c/f) - c = -fract(c/f)
1456  *
1457  * and we can add the same value to the original constraint.
1458  *
1459  * In the example, the compressed space only contains "j",
1460  * and the inequality translates to
1461  *
1462  *              3 * j - 1 >= 0
1463  *
1464  * We add -fract(-1/3) = -2 to the original constraint to obtain
1465  *
1466  *              i - 3 >= 0
1467  */
1468 static struct isl_basic_set *normalize_constraints_in_compressed_space(
1469         struct isl_basic_set *bset)
1470 {
1471         int i;
1472         unsigned total;
1473         struct isl_mat *B, *C;
1474         isl_int gcd;
1475
1476         if (!bset)
1477                 return NULL;
1478
1479         if (ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL))
1480                 return bset;
1481
1482         if (!bset->n_ineq)
1483                 return bset;
1484
1485         bset = isl_basic_set_cow(bset);
1486         if (!bset)
1487                 return NULL;
1488
1489         total = isl_basic_set_total_dim(bset);
1490         B = isl_mat_sub_alloc(bset->ctx, bset->eq, 0, bset->n_eq, 0, 1 + total);
1491         C = isl_mat_variable_compression(B, NULL);
1492         if (!C)
1493                 return bset;
1494         if (C->n_col == 0) {
1495                 isl_mat_free(C);
1496                 return isl_basic_set_set_to_empty(bset);
1497         }
1498         B = isl_mat_sub_alloc(bset->ctx, bset->ineq,
1499                                                 0, bset->n_ineq, 0, 1 + total);
1500         C = isl_mat_product(B, C);
1501         if (!C)
1502                 return bset;
1503
1504         isl_int_init(gcd);
1505         for (i = 0; i < bset->n_ineq; ++i) {
1506                 isl_seq_gcd(C->row[i] + 1, C->n_col - 1, &gcd);
1507                 if (isl_int_is_one(gcd))
1508                         continue;
1509                 isl_int_fdiv_r(C->row[i][0], C->row[i][0], gcd);
1510                 isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], C->row[i][0]);
1511         }
1512         isl_int_clear(gcd);
1513
1514         isl_mat_free(C);
1515
1516         return bset;
1517 }
1518
1519 /* Remove all information from bset that is redundant in the context
1520  * of context.  In particular, equalities that are linear combinations
1521  * of those in context are removed.  Then the inequalities that are
1522  * redundant in the context of the equalities and inequalities of
1523  * context are removed.
1524  *
1525  * We first simplify the constraints of "bset" in the context of the
1526  * equalities of "context".
1527  * Then we simplify the inequalities of the context in the context
1528  * of the equalities of bset and remove the inequalities from "bset"
1529  * that are obviously redundant with respect to some inequality in "context".
1530  *
1531  * If there are any inequalities left, we construct a tableau for
1532  * the context and then add the inequalities of "bset".
1533  * Before adding these equalities, we freeze all constraints such that
1534  * they won't be considered redundant in terms of the constraints of "bset".
1535  * Then we detect all equalities and redundant constraints (among the
1536  * constraints that weren't frozen) and update bset according to the results.
1537  * We have to be careful here because we don't want any of the context
1538  * constraints to remain and because we haven't added the equalities of "bset"
1539  * to the tableau so we temporarily have to pretend that there were no
1540  * equalities.
1541  */
1542 static struct isl_basic_set *uset_gist(struct isl_basic_set *bset,
1543         struct isl_basic_set *context)
1544 {
1545         int i;
1546         struct isl_tab *tab;
1547         unsigned context_ineq;
1548         struct isl_basic_set *combined = NULL;
1549
1550         if (!context || !bset)
1551                 goto error;
1552
1553         if (context->n_eq > 0)
1554                 bset = isl_basic_set_reduce_using_equalities(bset,
1555                                         isl_basic_set_copy(context));
1556         if (!bset)
1557                 goto error;
1558         if (isl_basic_set_fast_is_empty(bset))
1559                 goto done;
1560         if (!bset->n_ineq)
1561                 goto done;
1562
1563         if (bset->n_eq > 0) {
1564                 struct isl_basic_set *affine_hull;
1565                 affine_hull = isl_basic_set_copy(bset);
1566                 affine_hull = isl_basic_set_cow(affine_hull);
1567                 if (!affine_hull)
1568                         goto error;
1569                 isl_basic_set_free_inequality(affine_hull, affine_hull->n_ineq);
1570                 context = isl_basic_set_intersect(context, affine_hull);
1571                 context = isl_basic_set_gauss(context, NULL);
1572                 context = normalize_constraints_in_compressed_space(context);
1573         }
1574         if (!context)
1575                 goto error;
1576         if (ISL_F_ISSET(context, ISL_BASIC_SET_EMPTY)) {
1577                 isl_basic_set_free(bset);
1578                 return context;
1579         }
1580         if (!context->n_ineq)
1581                 goto done;
1582         bset = remove_shifted_constraints(bset, context);
1583         if (!bset->n_ineq)
1584                 goto done;
1585         context_ineq = context->n_ineq;
1586         combined = isl_basic_set_cow(isl_basic_set_copy(context));
1587         if (isl_basic_set_free_equality(combined, context->n_eq) < 0)
1588                 goto error;
1589         combined = isl_basic_set_extend_constraints(combined,
1590                                                     bset->n_eq, bset->n_ineq);
1591         tab = isl_tab_from_basic_set(combined);
1592         if (!tab)
1593                 goto error;
1594         for (i = 0; i < context_ineq; ++i)
1595                 if (isl_tab_freeze_constraint(tab, i) < 0)
1596                         goto error;
1597         tab = isl_tab_extend(tab, bset->n_ineq);
1598         if (!tab)
1599                 goto error;
1600         for (i = 0; i < bset->n_ineq; ++i)
1601                 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
1602                         goto error;
1603         bset = isl_basic_set_add_constraints(combined, bset, 0);
1604         tab = isl_tab_detect_implicit_equalities(tab);
1605         if (isl_tab_detect_redundant(tab) < 0) {
1606                 isl_tab_free(tab);
1607                 goto error2;
1608         }
1609         for (i = 0; i < context_ineq; ++i) {
1610                 tab->con[i].is_zero = 0;
1611                 tab->con[i].is_redundant = 1;
1612         }
1613         bset = isl_basic_set_update_from_tab(bset, tab);
1614         isl_tab_free(tab);
1615         ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
1616         ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
1617 done:
1618         bset = isl_basic_set_simplify(bset);
1619         bset = isl_basic_set_finalize(bset);
1620         isl_basic_set_free(context);
1621         return bset;
1622 error:
1623         isl_basic_set_free(combined);
1624 error2:
1625         isl_basic_set_free(bset);
1626         isl_basic_set_free(context);
1627         return NULL;
1628 }
1629
1630 /* Normalize the divs in "bmap" in the context of the equalities in "context".
1631  * We simply add the equalities in context to bmap and then do a regular
1632  * div normalizations.  Better results can be obtained by normalizing
1633  * only the divs in bmap than do not also appear in context.
1634  * We need to be careful to reduce the divs using the equalities
1635  * so that later calls to isl_basic_map_overlying_set wouldn't introduce
1636  * spurious constraints.
1637  */
1638 static struct isl_basic_map *normalize_divs_in_context(
1639         struct isl_basic_map *bmap, struct isl_basic_map *context)
1640 {
1641         int i;
1642         unsigned total_context;
1643         int div_eq;
1644
1645         div_eq = n_pure_div_eq(bmap);
1646         if (div_eq == 0)
1647                 return bmap;
1648
1649         if (context->n_div > 0)
1650                 bmap = isl_basic_map_align_divs(bmap, context);
1651
1652         total_context = isl_basic_map_total_dim(context);
1653         bmap = isl_basic_map_extend_constraints(bmap, context->n_eq, 0);
1654         for (i = 0; i < context->n_eq; ++i) {
1655                 int k;
1656                 k = isl_basic_map_alloc_equality(bmap);
1657                 isl_seq_cpy(bmap->eq[k], context->eq[i], 1 + total_context);
1658                 isl_seq_clr(bmap->eq[k] + 1 + total_context,
1659                                 isl_basic_map_total_dim(bmap) - total_context);
1660         }
1661         bmap = isl_basic_map_gauss(bmap, NULL);
1662         bmap = normalize_divs(bmap, NULL);
1663         bmap = isl_basic_map_gauss(bmap, NULL);
1664         return bmap;
1665 }
1666
1667 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
1668         struct isl_basic_map *context)
1669 {
1670         struct isl_basic_set *bset;
1671
1672         if (!bmap || !context)
1673                 goto error;
1674
1675         if (isl_basic_map_is_universe(context)) {
1676                 isl_basic_map_free(context);
1677                 return bmap;
1678         }
1679         if (isl_basic_map_is_universe(bmap)) {
1680                 isl_basic_map_free(context);
1681                 return bmap;
1682         }
1683         if (isl_basic_map_fast_is_empty(context)) {
1684                 struct isl_dim *dim = isl_dim_copy(bmap->dim);
1685                 isl_basic_map_free(context);
1686                 isl_basic_map_free(bmap);
1687                 return isl_basic_map_universe(dim);
1688         }
1689         if (isl_basic_map_fast_is_empty(bmap)) {
1690                 isl_basic_map_free(context);
1691                 return bmap;
1692         }
1693
1694         bmap = isl_basic_map_convex_hull(bmap);
1695         context = isl_basic_map_convex_hull(context);
1696
1697         if (context->n_eq)
1698                 bmap = normalize_divs_in_context(bmap, context);
1699
1700         context = isl_basic_map_align_divs(context, bmap);
1701         bmap = isl_basic_map_align_divs(bmap, context);
1702
1703         bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
1704                          isl_basic_map_underlying_set(context));
1705
1706         return isl_basic_map_overlying_set(bset, bmap);
1707 error:
1708         isl_basic_map_free(bmap);
1709         isl_basic_map_free(context);
1710         return NULL;
1711 }
1712
1713 /*
1714  * Assumes context has no implicit divs.
1715  */
1716 struct isl_map *isl_map_gist(struct isl_map *map, struct isl_basic_map *context)
1717 {
1718         int i;
1719
1720         if (!map || !context)
1721                 goto error;;
1722
1723         if (isl_basic_map_is_universe(context)) {
1724                 isl_basic_map_free(context);
1725                 return map;
1726         }
1727         if (isl_basic_map_fast_is_empty(context)) {
1728                 struct isl_dim *dim = isl_dim_copy(map->dim);
1729                 isl_basic_map_free(context);
1730                 isl_map_free(map);
1731                 return isl_map_universe(dim);
1732         }
1733
1734         context = isl_basic_map_convex_hull(context);
1735         map = isl_map_cow(map);
1736         if (!map || !context)
1737                 goto error;;
1738         isl_assert(map->ctx, isl_dim_equal(map->dim, context->dim), goto error);
1739         map = isl_map_compute_divs(map);
1740         for (i = 0; i < map->n; ++i)
1741                 context = isl_basic_map_align_divs(context, map->p[i]);
1742         for (i = 0; i < map->n; ++i) {
1743                 map->p[i] = isl_basic_map_gist(map->p[i],
1744                                                 isl_basic_map_copy(context));
1745                 if (!map->p[i])
1746                         goto error;
1747         }
1748         isl_basic_map_free(context);
1749         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
1750         return map;
1751 error:
1752         isl_map_free(map);
1753         isl_basic_map_free(context);
1754         return NULL;
1755 }
1756
1757 struct isl_basic_set *isl_basic_set_gist(struct isl_basic_set *bset,
1758                                                 struct isl_basic_set *context)
1759 {
1760         return (struct isl_basic_set *)isl_basic_map_gist(
1761                 (struct isl_basic_map *)bset, (struct isl_basic_map *)context);
1762 }
1763
1764 struct isl_set *isl_set_gist(struct isl_set *set, struct isl_basic_set *context)
1765 {
1766         return (struct isl_set *)isl_map_gist((struct isl_map *)set,
1767                                         (struct isl_basic_map *)context);
1768 }
1769
1770 /* Quick check to see if two basic maps are disjoint.
1771  * In particular, we reduce the equalities and inequalities of
1772  * one basic map in the context of the equalities of the other
1773  * basic map and check if we get a contradiction.
1774  */
1775 int isl_basic_map_fast_is_disjoint(struct isl_basic_map *bmap1,
1776         struct isl_basic_map *bmap2)
1777 {
1778         struct isl_vec *v = NULL;
1779         int *elim = NULL;
1780         unsigned total;
1781         int i;
1782
1783         if (!bmap1 || !bmap2)
1784                 return -1;
1785         isl_assert(bmap1->ctx, isl_dim_equal(bmap1->dim, bmap2->dim),
1786                         return -1);
1787         if (bmap1->n_div || bmap2->n_div)
1788                 return 0;
1789         if (!bmap1->n_eq && !bmap2->n_eq)
1790                 return 0;
1791
1792         total = isl_dim_total(bmap1->dim);
1793         if (total == 0)
1794                 return 0;
1795         v = isl_vec_alloc(bmap1->ctx, 1 + total);
1796         if (!v)
1797                 goto error;
1798         elim = isl_alloc_array(bmap1->ctx, int, total);
1799         if (!elim)
1800                 goto error;
1801         compute_elimination_index(bmap1, elim);
1802         for (i = 0; i < bmap2->n_eq; ++i) {
1803                 int reduced;
1804                 reduced = reduced_using_equalities(v->block.data, bmap2->eq[i],
1805                                                         bmap1, elim);
1806                 if (reduced && !isl_int_is_zero(v->block.data[0]) &&
1807                     isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1808                         goto disjoint;
1809         }
1810         for (i = 0; i < bmap2->n_ineq; ++i) {
1811                 int reduced;
1812                 reduced = reduced_using_equalities(v->block.data,
1813                                                 bmap2->ineq[i], bmap1, elim);
1814                 if (reduced && isl_int_is_neg(v->block.data[0]) &&
1815                     isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1816                         goto disjoint;
1817         }
1818         compute_elimination_index(bmap2, elim);
1819         for (i = 0; i < bmap1->n_ineq; ++i) {
1820                 int reduced;
1821                 reduced = reduced_using_equalities(v->block.data,
1822                                                 bmap1->ineq[i], bmap2, elim);
1823                 if (reduced && isl_int_is_neg(v->block.data[0]) &&
1824                     isl_seq_first_non_zero(v->block.data + 1, total) == -1)
1825                         goto disjoint;
1826         }
1827         isl_vec_free(v);
1828         free(elim);
1829         return 0;
1830 disjoint:
1831         isl_vec_free(v);
1832         free(elim);
1833         return 1;
1834 error:
1835         isl_vec_free(v);
1836         free(elim);
1837         return -1;
1838 }
1839
1840 int isl_basic_set_fast_is_disjoint(struct isl_basic_set *bset1,
1841         struct isl_basic_set *bset2)
1842 {
1843         return isl_basic_map_fast_is_disjoint((struct isl_basic_map *)bset1,
1844                                               (struct isl_basic_map *)bset2);
1845 }
1846
1847 int isl_map_fast_is_disjoint(struct isl_map *map1, struct isl_map *map2)
1848 {
1849         int i, j;
1850
1851         if (!map1 || !map2)
1852                 return -1;
1853
1854         if (isl_map_fast_is_equal(map1, map2))
1855                 return 0;
1856
1857         for (i = 0; i < map1->n; ++i) {
1858                 for (j = 0; j < map2->n; ++j) {
1859                         int d = isl_basic_map_fast_is_disjoint(map1->p[i],
1860                                                                map2->p[j]);
1861                         if (d != 1)
1862                                 return d;
1863                 }
1864         }
1865         return 1;
1866 }
1867
1868 int isl_set_fast_is_disjoint(struct isl_set *set1, struct isl_set *set2)
1869 {
1870         return isl_map_fast_is_disjoint((struct isl_map *)set1,
1871                                         (struct isl_map *)set2);
1872 }
1873
1874 /* Check if we can combine a given div with lower bound l and upper
1875  * bound u with some other div and if so return that other div.
1876  * Otherwise return -1.
1877  *
1878  * We first check that
1879  *      - the bounds are opposites of each other (except for the constant
1880  *        term)
1881  *      - the bounds do not reference any other div
1882  *      - no div is defined in terms of this div
1883  *
1884  * Let m be the size of the range allowed on the div by the bounds.
1885  * That is, the bounds are of the form
1886  *
1887  *      e <= a <= e + m - 1
1888  *
1889  * with e some expression in the other variables.
1890  * We look for another div b such that no third div is defined in terms
1891  * of this second div b and such that in any constraint that contains
1892  * a (except for the given lower and upper bound), also contains b
1893  * with a coefficient that is m times that of b.
1894  * That is, all constraints (execpt for the lower and upper bound)
1895  * are of the form
1896  *
1897  *      e + f (a + m b) >= 0
1898  *
1899  * If so, we return b so that "a + m b" can be replaced by
1900  * a single div "c = a + m b".
1901  */
1902 static int div_find_coalesce(struct isl_basic_map *bmap, int *pairs,
1903         unsigned div, unsigned l, unsigned u)
1904 {
1905         int i, j;
1906         unsigned dim;
1907         int coalesce = -1;
1908
1909         if (bmap->n_div <= 1)
1910                 return -1;
1911         dim = isl_dim_total(bmap->dim);
1912         if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim, div) != -1)
1913                 return -1;
1914         if (isl_seq_first_non_zero(bmap->ineq[l] + 1 + dim + div + 1,
1915                                    bmap->n_div - div - 1) != -1)
1916                 return -1;
1917         if (!isl_seq_is_neg(bmap->ineq[l] + 1, bmap->ineq[u] + 1,
1918                             dim + bmap->n_div))
1919                 return -1;
1920
1921         for (i = 0; i < bmap->n_div; ++i) {
1922                 if (isl_int_is_zero(bmap->div[i][0]))
1923                         continue;
1924                 if (!isl_int_is_zero(bmap->div[i][1 + 1 + dim + div]))
1925                         return -1;
1926         }
1927
1928         isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
1929         if (isl_int_is_neg(bmap->ineq[l][0])) {
1930                 isl_int_sub(bmap->ineq[l][0],
1931                             bmap->ineq[l][0], bmap->ineq[u][0]);
1932                 bmap = isl_basic_map_copy(bmap);
1933                 bmap = isl_basic_map_set_to_empty(bmap);
1934                 isl_basic_map_free(bmap);
1935                 return -1;
1936         }
1937         isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
1938         for (i = 0; i < bmap->n_div; ++i) {
1939                 if (i == div)
1940                         continue;
1941                 if (!pairs[i])
1942                         continue;
1943                 for (j = 0; j < bmap->n_div; ++j) {
1944                         if (isl_int_is_zero(bmap->div[j][0]))
1945                                 continue;
1946                         if (!isl_int_is_zero(bmap->div[j][1 + 1 + dim + i]))
1947                                 break;
1948                 }
1949                 if (j < bmap->n_div)
1950                         continue;
1951                 for (j = 0; j < bmap->n_ineq; ++j) {
1952                         int valid;
1953                         if (j == l || j == u)
1954                                 continue;
1955                         if (isl_int_is_zero(bmap->ineq[j][1 + dim + div]))
1956                                 continue;
1957                         if (isl_int_is_zero(bmap->ineq[j][1 + dim + i]))
1958                                 break;
1959                         isl_int_mul(bmap->ineq[j][1 + dim + div],
1960                                     bmap->ineq[j][1 + dim + div],
1961                                     bmap->ineq[l][0]);
1962                         valid = isl_int_eq(bmap->ineq[j][1 + dim + div],
1963                                            bmap->ineq[j][1 + dim + i]);
1964                         isl_int_divexact(bmap->ineq[j][1 + dim + div],
1965                                          bmap->ineq[j][1 + dim + div],
1966                                          bmap->ineq[l][0]);
1967                         if (!valid)
1968                                 break;
1969                 }
1970                 if (j < bmap->n_ineq)
1971                         continue;
1972                 coalesce = i;
1973                 break;
1974         }
1975         isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
1976         isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
1977         return coalesce;
1978 }
1979
1980 /* Given a lower and an upper bound on div i, construct an inequality
1981  * that when nonnegative ensures that this pair of bounds always allows
1982  * for an integer value of the given div.
1983  * The lower bound is inequality l, while the upper bound is inequality u.
1984  * The constructed inequality is stored in ineq.
1985  * g, fl, fu are temporary scalars.
1986  *
1987  * Let the upper bound be
1988  *
1989  *      -n_u a + e_u >= 0
1990  *
1991  * and the lower bound
1992  *
1993  *      n_l a + e_l >= 0
1994  *
1995  * Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
1996  * We have
1997  *
1998  *      - f_u e_l <= f_u f_l g a <= f_l e_u
1999  *
2000  * Since all variables are integer valued, this is equivalent to
2001  *
2002  *      - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
2003  *
2004  * If this interval is at least f_u f_l g, then it contains at least
2005  * one integer value for a.
2006  * That is, the test constraint is
2007  *
2008  *      f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
2009  */
2010 static void construct_test_ineq(struct isl_basic_map *bmap, int i,
2011         int l, int u, isl_int *ineq, isl_int g, isl_int fl, isl_int fu)
2012 {
2013         unsigned dim;
2014         dim = isl_dim_total(bmap->dim);
2015
2016         isl_int_gcd(g, bmap->ineq[l][1 + dim + i], bmap->ineq[u][1 + dim + i]);
2017         isl_int_divexact(fl, bmap->ineq[l][1 + dim + i], g);
2018         isl_int_divexact(fu, bmap->ineq[u][1 + dim + i], g);
2019         isl_int_neg(fu, fu);
2020         isl_seq_combine(ineq, fl, bmap->ineq[u], fu, bmap->ineq[l],
2021                         1 + dim + bmap->n_div);
2022         isl_int_add(ineq[0], ineq[0], fl);
2023         isl_int_add(ineq[0], ineq[0], fu);
2024         isl_int_sub_ui(ineq[0], ineq[0], 1);
2025         isl_int_mul(g, g, fl);
2026         isl_int_mul(g, g, fu);
2027         isl_int_sub(ineq[0], ineq[0], g);
2028 }
2029
2030 /* Remove more kinds of divs that are not strictly needed.
2031  * In particular, if all pairs of lower and upper bounds on a div
2032  * are such that they allow at least one integer value of the div,
2033  * the we can eliminate the div using Fourier-Motzkin without
2034  * introducing any spurious solutions.
2035  */
2036 static struct isl_basic_map *drop_more_redundant_divs(
2037         struct isl_basic_map *bmap, int *pairs, int n)
2038 {
2039         struct isl_tab *tab = NULL;
2040         struct isl_vec *vec = NULL;
2041         unsigned dim;
2042         int remove = -1;
2043         isl_int g, fl, fu;
2044
2045         isl_int_init(g);
2046         isl_int_init(fl);
2047         isl_int_init(fu);
2048
2049         if (!bmap)
2050                 goto error;
2051
2052         dim = isl_dim_total(bmap->dim);
2053         vec = isl_vec_alloc(bmap->ctx, 1 + dim + bmap->n_div);
2054         if (!vec)
2055                 goto error;
2056
2057         tab = isl_tab_from_basic_map(bmap);
2058
2059         while (n > 0) {
2060                 int i, l, u;
2061                 int best = -1;
2062                 enum isl_lp_result res;
2063
2064                 for (i = 0; i < bmap->n_div; ++i) {
2065                         if (!pairs[i])
2066                                 continue;
2067                         if (best >= 0 && pairs[best] <= pairs[i])
2068                                 continue;
2069                         best = i;
2070                 }
2071
2072                 i = best;
2073                 for (l = 0; l < bmap->n_ineq; ++l) {
2074                         if (!isl_int_is_pos(bmap->ineq[l][1 + dim + i]))
2075                                 continue;
2076                         for (u = 0; u < bmap->n_ineq; ++u) {
2077                                 if (!isl_int_is_neg(bmap->ineq[u][1 + dim + i]))
2078                                         continue;
2079                                 construct_test_ineq(bmap, i, l, u,
2080                                                     vec->el, g, fl, fu);
2081                                 res = isl_tab_min(tab, vec->el,
2082                                                   bmap->ctx->one, &g, NULL, 0);
2083                                 if (res == isl_lp_error)
2084                                         goto error;
2085                                 if (res == isl_lp_empty) {
2086                                         bmap = isl_basic_map_set_to_empty(bmap);
2087                                         break;
2088                                 }
2089                                 if (res != isl_lp_ok || isl_int_is_neg(g))
2090                                         break;
2091                         }
2092                         if (u < bmap->n_ineq)
2093                                 break;
2094                 }
2095                 if (l == bmap->n_ineq) {
2096                         remove = i;
2097                         break;
2098                 }
2099                 pairs[i] = 0;
2100                 --n;
2101         }
2102
2103         isl_tab_free(tab);
2104         isl_vec_free(vec);
2105
2106         isl_int_clear(g);
2107         isl_int_clear(fl);
2108         isl_int_clear(fu);
2109
2110         free(pairs);
2111
2112         if (remove < 0)
2113                 return bmap;
2114
2115         bmap = isl_basic_map_remove(bmap, isl_dim_div, remove, 1);
2116         return isl_basic_map_drop_redundant_divs(bmap);
2117 error:
2118         free(pairs);
2119         isl_basic_map_free(bmap);
2120         isl_tab_free(tab);
2121         isl_vec_free(vec);
2122         isl_int_clear(g);
2123         isl_int_clear(fl);
2124         isl_int_clear(fu);
2125         return NULL;
2126 }
2127
2128 /* Given a pair of divs div1 and div2 such that, expect for the lower bound l
2129  * and the upper bound u, div1 always occurs together with div2 in the form
2130  * (div1 + m div2), where m is the constant range on the variable div1
2131  * allowed by l and u, replace the pair div1 and div2 by a single
2132  * div that is equal to div1 + m div2.
2133  *
2134  * The new div will appear in the location that contains div2.
2135  * We need to modify all constraints that contain
2136  * div2 = (div - div1) / m
2137  * (If a constraint does not contain div2, it will also not contain div1.)
2138  * If the constraint also contains div1, then we know they appear
2139  * as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
2140  * i.e., the coefficient of div is f.
2141  *
2142  * Otherwise, we first need to introduce div1 into the constraint.
2143  * Let the l be
2144  *
2145  *      div1 + f >=0
2146  *
2147  * and u
2148  *
2149  *      -div1 + f' >= 0
2150  *
2151  * A lower bound on div2
2152  *
2153  *      n div2 + t >= 0
2154  *
2155  * can be replaced by
2156  *
2157  *      (n * (m div 2 + div1) + m t + n f)/g >= 0
2158  *
2159  * with g = gcd(m,n).
2160  * An upper bound
2161  *
2162  *      -n div2 + t >= 0
2163  *
2164  * can be replaced by
2165  *
2166  *      (-n * (m div2 + div1) + m t + n f')/g >= 0
2167  *
2168  * These constraint are those that we would obtain from eliminating
2169  * div1 using Fourier-Motzkin.
2170  *
2171  * After all constraints have been modified, we drop the lower and upper
2172  * bound and then drop div1.
2173  */
2174 static struct isl_basic_map *coalesce_divs(struct isl_basic_map *bmap,
2175         unsigned div1, unsigned div2, unsigned l, unsigned u)
2176 {
2177         isl_int a;
2178         isl_int b;
2179         isl_int m;
2180         unsigned dim, total;
2181         int i;
2182
2183         dim = isl_dim_total(bmap->dim);
2184         total = 1 + dim + bmap->n_div;
2185
2186         isl_int_init(a);
2187         isl_int_init(b);
2188         isl_int_init(m);
2189         isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
2190         isl_int_add_ui(m, m, 1);
2191
2192         for (i = 0; i < bmap->n_ineq; ++i) {
2193                 if (i == l || i == u)
2194                         continue;
2195                 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div2]))
2196                         continue;
2197                 if (isl_int_is_zero(bmap->ineq[i][1 + dim + div1])) {
2198                         isl_int_gcd(b, m, bmap->ineq[i][1 + dim + div2]);
2199                         isl_int_divexact(a, m, b);
2200                         isl_int_divexact(b, bmap->ineq[i][1 + dim + div2], b);
2201                         if (isl_int_is_pos(b)) {
2202                                 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
2203                                                 b, bmap->ineq[l], total);
2204                         } else {
2205                                 isl_int_neg(b, b);
2206                                 isl_seq_combine(bmap->ineq[i], a, bmap->ineq[i],
2207                                                 b, bmap->ineq[u], total);
2208                         }
2209                 }
2210                 isl_int_set(bmap->ineq[i][1 + dim + div2],
2211                             bmap->ineq[i][1 + dim + div1]);
2212                 isl_int_set_si(bmap->ineq[i][1 + dim + div1], 0);
2213         }
2214
2215         isl_int_clear(a);
2216         isl_int_clear(b);
2217         isl_int_clear(m);
2218         if (l > u) {
2219                 isl_basic_map_drop_inequality(bmap, l);
2220                 isl_basic_map_drop_inequality(bmap, u);
2221         } else {
2222                 isl_basic_map_drop_inequality(bmap, u);
2223                 isl_basic_map_drop_inequality(bmap, l);
2224         }
2225         bmap = isl_basic_map_drop_div(bmap, div1);
2226         return bmap;
2227 }
2228
2229 /* First check if we can coalesce any pair of divs and
2230  * then continue with dropping more redundant divs.
2231  *
2232  * We loop over all pairs of lower and upper bounds on a div
2233  * with coefficient 1 and -1, respectively, check if there
2234  * is any other div "c" with which we can coalesce the div
2235  * and if so, perform the coalescing.
2236  */
2237 static struct isl_basic_map *coalesce_or_drop_more_redundant_divs(
2238         struct isl_basic_map *bmap, int *pairs, int n)
2239 {
2240         int i, l, u;
2241         unsigned dim;
2242
2243         dim = isl_dim_total(bmap->dim);
2244
2245         for (i = 0; i < bmap->n_div; ++i) {
2246                 if (!pairs[i])
2247                         continue;
2248                 for (l = 0; l < bmap->n_ineq; ++l) {
2249                         if (!isl_int_is_one(bmap->ineq[l][1 + dim + i]))
2250                                 continue;
2251                         for (u = 0; u < bmap->n_ineq; ++u) {
2252                                 int c;
2253
2254                                 if (!isl_int_is_negone(bmap->ineq[u][1+dim+i]))
2255                                         continue;
2256                                 c = div_find_coalesce(bmap, pairs, i, l, u);
2257                                 if (c < 0)
2258                                         continue;
2259                                 free(pairs);
2260                                 bmap = coalesce_divs(bmap, i, c, l, u);
2261                                 return isl_basic_map_drop_redundant_divs(bmap);
2262                         }
2263                 }
2264         }
2265
2266         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
2267                 return bmap;
2268
2269         return drop_more_redundant_divs(bmap, pairs, n);
2270 }
2271
2272 /* Remove divs that are not strictly needed.
2273  * In particular, if a div only occurs positively (or negatively)
2274  * in constraints, then it can simply be dropped.
2275  * Also, if a div occurs only occurs in two constraints and if moreover
2276  * those two constraints are opposite to each other, except for the constant
2277  * term and if the sum of the constant terms is such that for any value
2278  * of the other values, there is always at least one integer value of the
2279  * div, i.e., if one plus this sum is greater than or equal to
2280  * the (absolute value) of the coefficent of the div in the constraints,
2281  * then we can also simply drop the div.
2282  *
2283  * If any divs are left after these simple checks then we move on
2284  * to more complicated cases in drop_more_redundant_divs.
2285  */
2286 struct isl_basic_map *isl_basic_map_drop_redundant_divs(
2287         struct isl_basic_map *bmap)
2288 {
2289         int i, j;
2290         unsigned off;
2291         int *pairs = NULL;
2292         int n = 0;
2293
2294         if (!bmap)
2295                 goto error;
2296
2297         off = isl_dim_total(bmap->dim);
2298         pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
2299         if (!pairs)
2300                 goto error;
2301
2302         for (i = 0; i < bmap->n_div; ++i) {
2303                 int pos, neg;
2304                 int last_pos, last_neg;
2305                 int redundant;
2306                 int defined;
2307
2308                 defined = !isl_int_is_zero(bmap->div[i][0]);
2309                 for (j = 0; j < bmap->n_eq; ++j)
2310                         if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
2311                                 break;
2312                 if (j < bmap->n_eq)
2313                         continue;
2314                 ++n;
2315                 pos = neg = 0;
2316                 for (j = 0; j < bmap->n_ineq; ++j) {
2317                         if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
2318                                 last_pos = j;
2319                                 ++pos;
2320                         }
2321                         if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
2322                                 last_neg = j;
2323                                 ++neg;
2324                         }
2325                 }
2326                 pairs[i] = pos * neg;
2327                 if (pairs[i] == 0) {
2328                         for (j = bmap->n_ineq - 1; j >= 0; --j)
2329                                 if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
2330                                         isl_basic_map_drop_inequality(bmap, j);
2331                         bmap = isl_basic_map_drop_div(bmap, i);
2332                         free(pairs);
2333                         return isl_basic_map_drop_redundant_divs(bmap);
2334                 }
2335                 if (pairs[i] != 1)
2336                         continue;
2337                 if (!isl_seq_is_neg(bmap->ineq[last_pos] + 1,
2338                                     bmap->ineq[last_neg] + 1,
2339                                     off + bmap->n_div))
2340                         continue;
2341
2342                 isl_int_add(bmap->ineq[last_pos][0],
2343                             bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
2344                 isl_int_add_ui(bmap->ineq[last_pos][0],
2345                                bmap->ineq[last_pos][0], 1);
2346                 redundant = isl_int_ge(bmap->ineq[last_pos][0],
2347                                 bmap->ineq[last_pos][1+off+i]);
2348                 isl_int_sub_ui(bmap->ineq[last_pos][0],
2349                                bmap->ineq[last_pos][0], 1);
2350                 isl_int_sub(bmap->ineq[last_pos][0],
2351                             bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
2352                 if (!redundant) {
2353                         if (defined ||
2354                             !ok_to_set_div_from_bound(bmap, i, last_pos)) {
2355                                 pairs[i] = 0;
2356                                 --n;
2357                                 continue;
2358                         }
2359                         bmap = set_div_from_lower_bound(bmap, i, last_pos);
2360                         bmap = isl_basic_map_simplify(bmap);
2361                         free(pairs);
2362                         return isl_basic_map_drop_redundant_divs(bmap);
2363                 }
2364                 if (last_pos > last_neg) {
2365                         isl_basic_map_drop_inequality(bmap, last_pos);
2366                         isl_basic_map_drop_inequality(bmap, last_neg);
2367                 } else {
2368                         isl_basic_map_drop_inequality(bmap, last_neg);
2369                         isl_basic_map_drop_inequality(bmap, last_pos);
2370                 }
2371                 bmap = isl_basic_map_drop_div(bmap, i);
2372                 free(pairs);
2373                 return isl_basic_map_drop_redundant_divs(bmap);
2374         }
2375
2376         if (n > 0)
2377                 return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
2378
2379         free(pairs);
2380         return bmap;
2381 error:
2382         free(pairs);
2383         isl_basic_map_free(bmap);
2384         return NULL;
2385 }
2386
2387 struct isl_basic_set *isl_basic_set_drop_redundant_divs(
2388         struct isl_basic_set *bset)
2389 {
2390         return (struct isl_basic_set *)
2391             isl_basic_map_drop_redundant_divs((struct isl_basic_map *)bset);
2392 }
2393
2394 struct isl_map *isl_map_drop_redundant_divs(struct isl_map *map)
2395 {
2396         int i;
2397
2398         if (!map)
2399                 return NULL;
2400         for (i = 0; i < map->n; ++i) {
2401                 map->p[i] = isl_basic_map_drop_redundant_divs(map->p[i]);
2402                 if (!map->p[i])
2403                         goto error;
2404         }
2405         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
2406         return map;
2407 error:
2408         isl_map_free(map);
2409         return NULL;
2410 }
2411
2412 struct isl_set *isl_set_drop_redundant_divs(struct isl_set *set)
2413 {
2414         return (struct isl_set *)
2415             isl_map_drop_redundant_divs((struct isl_map *)set);
2416 }