isl_affine_hull.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_ctx.h"
  11 #include "isl_seq.h"
  12 #include "isl_set.h"
  13 #include "isl_lp.h"
  14 #include "isl_map.h"
  15 #include "isl_map_private.h"
  16 #include "isl_equalities.h"
  17 #include "isl_sample.h"
  18 #include "isl_tab.h"
  19
  20 struct isl_basic_map *isl_basic_map_implicit_equalities(
  21                                                 struct isl_basic_map *bmap)
  22 {
  23         struct isl_tab *tab;
  24
  25         if (!bmap)
  26                 return bmap;
  27
  28         bmap = isl_basic_map_gauss(bmap, NULL);
  29         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
  30                 return bmap;
  31         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_IMPLICIT))
  32                 return bmap;
  33         if (bmap->n_ineq <= 1)
  34                 return bmap;
  35
  36         tab = isl_tab_from_basic_map(bmap);
  37         tab = isl_tab_detect_implicit_equalities(tab);
  38         bmap = isl_basic_map_update_from_tab(bmap, tab);
  39         isl_tab_free(tab);
  40         bmap = isl_basic_map_gauss(bmap, NULL);
  41         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
  42         return bmap;
  43 }
  44
  45 struct isl_basic_set *isl_basic_set_implicit_equalities(
  46                                                 struct isl_basic_set *bset)
  47 {
  48         return (struct isl_basic_set *)
  49                 isl_basic_map_implicit_equalities((struct isl_basic_map*)bset);
  50 }
  51
  52 struct isl_map *isl_map_implicit_equalities(struct isl_map *map)
  53 {
  54         int i;
  55
  56         if (!map)
  57                 return map;
  58
  59         for (i = 0; i < map->n; ++i) {
  60                 map->p[i] = isl_basic_map_implicit_equalities(map->p[i]);
  61                 if (!map->p[i])
  62                         goto error;
  63         }
  64
  65         return map;
  66 error:
  67         isl_map_free(map);
  68         return NULL;
  69 }
  70
  71 /* Make eq[row][col] of both bmaps equal so we can add the row
  72  * add the column to the common matrix.
  73  * Note that because of the echelon form, the columns of row row
  74  * after column col are zero.
  75  */
  76 static void set_common_multiple(
  77         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
  78         unsigned row, unsigned col)
  79 {
  80         isl_int m, c;
  81
  82         if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
  83                 return;
  84
  85         isl_int_init(c);
  86         isl_int_init(m);
  87         isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
  88         isl_int_divexact(c, m, bset1->eq[row][col]);
  89         isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
  90         isl_int_divexact(c, m, bset2->eq[row][col]);
  91         isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
  92         isl_int_clear(c);
  93         isl_int_clear(m);
  94 }
  95
  96 /* Delete a given equality, moving all the following equalities one up.
  97  */
  98 static void delete_row(struct isl_basic_set *bset, unsigned row)
  99 {
 100         isl_int *t;
 101         int r;
 102
 103         t = bset->eq[row];
 104         bset->n_eq--;
 105         for (r = row; r < bset->n_eq; ++r)
 106                 bset->eq[r] = bset->eq[r+1];
 107         bset->eq[bset->n_eq] = t;
 108 }
 109
 110 /* Make first row entries in column col of bset1 identical to
 111  * those of bset2, using the fact that entry bset1->eq[row][col]=a
 112  * is non-zero.  Initially, these elements of bset1 are all zero.
 113  * For each row i < row, we set
 114  *              A[i] = a * A[i] + B[i][col] * A[row]
 115  *              B[i] = a * B[i]
 116  * so that
 117  *              A[i][col] = B[i][col] = a * old(B[i][col])
 118  */
 119 static void construct_column(
 120         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 121         unsigned row, unsigned col)
 122 {
 123         int r;
 124         isl_int a;
 125         isl_int b;
 126         unsigned total;
 127
 128         isl_int_init(a);
 129         isl_int_init(b);
 130         total = 1 + isl_basic_set_n_dim(bset1);
 131         for (r = 0; r < row; ++r) {
 132                 if (isl_int_is_zero(bset2->eq[r][col]))
 133                         continue;
 134                 isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
 135                 isl_int_divexact(a, bset1->eq[row][col], b);
 136                 isl_int_divexact(b, bset2->eq[r][col], b);
 137                 isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
 138                                               b, bset1->eq[row], total);
 139                 isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 140         }
 141         isl_int_clear(a);
 142         isl_int_clear(b);
 143         delete_row(bset1, row);
 144 }
 145
 146 /* Make first row entries in column col of bset1 identical to
 147  * those of bset2, using only these entries of the two matrices.
 148  * Let t be the last row with different entries.
 149  * For each row i < t, we set
 150  *      A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
 151  *      B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
 152  * so that
 153  *      A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
 154  */
 155 static int transform_column(
 156         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 157         unsigned row, unsigned col)
 158 {
 159         int i, t;
 160         isl_int a, b, g;
 161         unsigned total;
 162
 163         for (t = row-1; t >= 0; --t)
 164                 if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 165                         break;
 166         if (t < 0)
 167                 return 0;
 168
 169         total = 1 + isl_basic_set_n_dim(bset1);
 170         isl_int_init(a);
 171         isl_int_init(b);
 172         isl_int_init(g);
 173         isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 174         for (i = 0; i < t; ++i) {
 175                 isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 176                 isl_int_gcd(g, a, b);
 177                 isl_int_divexact(a, a, g);
 178                 isl_int_divexact(g, b, g);
 179                 isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 180                                 total);
 181                 isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 182                                 total);
 183         }
 184         isl_int_clear(a);
 185         isl_int_clear(b);
 186         isl_int_clear(g);
 187         delete_row(bset1, t);
 188         delete_row(bset2, t);
 189         return 1;
 190 }
 191
 192 /* The implementation is based on Section 5.2 of Michael Karr,
 193  * "Affine Relationships Among Variables of a Program",
 194  * except that the echelon form we use starts from the last column
 195  * and that we are dealing with integer coefficients.
 196  */
 197 static struct isl_basic_set *affine_hull(
 198         struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 199 {
 200         unsigned total;
 201         int col;
 202         int row;
 203
 204         total = 1 + isl_basic_set_n_dim(bset1);
 205
 206         row = 0;
 207         for (col = total-1; col >= 0; --col) {
 208                 int is_zero1 = row >= bset1->n_eq ||
 209                         isl_int_is_zero(bset1->eq[row][col]);
 210                 int is_zero2 = row >= bset2->n_eq ||
 211                         isl_int_is_zero(bset2->eq[row][col]);
 212                 if (!is_zero1 && !is_zero2) {
 213                         set_common_multiple(bset1, bset2, row, col);
 214                         ++row;
 215                 } else if (!is_zero1 && is_zero2) {
 216                         construct_column(bset1, bset2, row, col);
 217                 } else if (is_zero1 && !is_zero2) {
 218                         construct_column(bset2, bset1, row, col);
 219                 } else {
 220                         if (transform_column(bset1, bset2, row, col))
 221                                 --row;
 222                 }
 223         }
 224         isl_basic_set_free(bset2);
 225         isl_assert(bset1->ctx, row == bset1->n_eq, goto error);
 226         bset1 = isl_basic_set_normalize_constraints(bset1);
 227         return bset1;
 228 error:
 229         isl_basic_set_free(bset1);
 230         return NULL;
 231 }
 232
 233 /* Find an integer point in the set represented by "tab"
 234  * that lies outside of the equality "eq" e(x) = 0.
 235  * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
 236  * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
 237  * The point, if found, is returned.
 238  * If no point can be found, a zero-length vector is returned.
 239  *
 240  * Before solving an ILP problem, we first check if simply
 241  * adding the normal of the constraint to one of the known
 242  * integer points in the basic set represented by "tab"
 243  * yields another point inside the basic set.
 244  *
 245  * The caller of this function ensures that the tableau is bounded or
 246  * that tab->basis and tab->n_unbounded have been set appropriately.
 247  */
 248 static struct isl_vec *outside_point(struct isl_tab *tab, isl_int *eq, int up)
 249 {
 250         struct isl_ctx *ctx;
 251         struct isl_vec *sample = NULL;
 252         struct isl_tab_undo *snap;
 253         unsigned dim;
 254
 255         if (!tab)
 256                 return NULL;
 257         ctx = tab->mat->ctx;
 258
 259         dim = tab->n_var;
 260         sample = isl_vec_alloc(ctx, 1 + dim);
 261         if (!sample)
 262                 return NULL;
 263         isl_int_set_si(sample->el[0], 1);
 264         isl_seq_combine(sample->el + 1,
 265                 ctx->one, tab->bmap->sample->el + 1,
 266                 up ? ctx->one : ctx->negone, eq + 1, dim);
 267         if (isl_basic_map_contains(tab->bmap, sample))
 268                 return sample;
 269         isl_vec_free(sample);
 270         sample = NULL;
 271
 272         snap = isl_tab_snap(tab);
 273
 274         if (!up)
 275                 isl_seq_neg(eq, eq, 1 + dim);
 276         isl_int_sub_ui(eq[0], eq[0], 1);
 277
 278         if (isl_tab_extend_cons(tab, 1) < 0)
 279                 goto error;
 280         if (isl_tab_add_ineq(tab, eq) < 0)
 281                 goto error;
 282
 283         sample = isl_tab_sample(tab);
 284
 285         isl_int_add_ui(eq[0], eq[0], 1);
 286         if (!up)
 287                 isl_seq_neg(eq, eq, 1 + dim);
 288
 289         if (isl_tab_rollback(tab, snap) < 0)
 290                 goto error;
 291
 292         return sample;
 293 error:
 294         isl_vec_free(sample);
 295         return NULL;
 296 }
 297
 298 struct isl_basic_set *isl_basic_set_recession_cone(struct isl_basic_set *bset)
 299 {
 300         int i;
 301
 302         bset = isl_basic_set_cow(bset);
 303         if (!bset)
 304                 return NULL;
 305         isl_assert(bset->ctx, bset->n_div == 0, goto error);
 306
 307         for (i = 0; i < bset->n_eq; ++i)
 308                 isl_int_set_si(bset->eq[i][0], 0);
 309
 310         for (i = 0; i < bset->n_ineq; ++i)
 311                 isl_int_set_si(bset->ineq[i][0], 0);
 312
 313         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
 314         return isl_basic_set_implicit_equalities(bset);
 315 error:
 316         isl_basic_set_free(bset);
 317         return NULL;
 318 }
 319
 320 __isl_give isl_set *isl_set_recession_cone(__isl_take isl_set *set)
 321 {
 322         int i;
 323
 324         if (!set)
 325                 return NULL;
 326         if (set->n == 0)
 327                 return set;
 328
 329         set = isl_set_remove_divs(set);
 330         set = isl_set_cow(set);
 331         if (!set)
 332                 return NULL;
 333
 334         for (i = 0; i < set->n; ++i) {
 335                 set->p[i] = isl_basic_set_recession_cone(set->p[i]);
 336                 if (!set->p[i])
 337                         goto error;
 338         }
 339
 340         return set;
 341 error:
 342         isl_set_free(set);
 343         return NULL;
 344 }
 345
 346 /* Extend an initial (under-)approximation of the affine hull of basic
 347  * set represented by the tableau "tab"
 348  * by looking for points that do not satisfy one of the equalities
 349  * in the current approximation and adding them to that approximation
 350  * until no such points can be found any more.
 351  *
 352  * The caller of this function ensures that "tab" is bounded or
 353  * that tab->basis and tab->n_unbounded have been set appropriately.
 354  */
 355 static struct isl_basic_set *extend_affine_hull(struct isl_tab *tab,
 356         struct isl_basic_set *hull)
 357 {
 358         int i, j;
 359         unsigned dim;
 360
 361         if (!tab || !hull)
 362                 goto error;
 363
 364         dim = tab->n_var;
 365
 366         if (isl_tab_extend_cons(tab, 2 * dim + 1) < 0)
 367                 goto error;
 368
 369         for (i = 0; i < dim; ++i) {
 370                 struct isl_vec *sample;
 371                 struct isl_basic_set *point;
 372                 for (j = 0; j < hull->n_eq; ++j) {
 373                         sample = outside_point(tab, hull->eq[j], 1);
 374                         if (!sample)
 375                                 goto error;
 376                         if (sample->size > 0)
 377                                 break;
 378                         isl_vec_free(sample);
 379                         sample = outside_point(tab, hull->eq[j], 0);
 380                         if (!sample)
 381                                 goto error;
 382                         if (sample->size > 0)
 383                                 break;
 384                         isl_vec_free(sample);
 385
 386                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 387                         if (!tab)
 388                                 goto error;
 389                 }
 390                 if (j == hull->n_eq)
 391                         break;
 392                 if (tab->samples)
 393                         tab = isl_tab_add_sample(tab, isl_vec_copy(sample));
 394                 if (!tab)
 395                         goto error;
 396                 point = isl_basic_set_from_vec(sample);
 397                 hull = affine_hull(hull, point);
 398         }
 399
 400         return hull;
 401 error:
 402         isl_basic_set_free(hull);
 403         return NULL;
 404 }
 405
 406 /* Drop all constraints in bset that involve any of the dimensions
 407  * first to first+n-1.
 408  */
 409 static struct isl_basic_set *drop_constraints_involving
 410         (struct isl_basic_set *bset, unsigned first, unsigned n)
 411 {
 412         int i;
 413
 414         if (!bset)
 415                 return NULL;
 416
 417         bset = isl_basic_set_cow(bset);
 418
 419         for (i = bset->n_eq - 1; i >= 0; --i) {
 420                 if (isl_seq_first_non_zero(bset->eq[i] + 1 + first, n) == -1)
 421                         continue;
 422                 isl_basic_set_drop_equality(bset, i);
 423         }
 424
 425         for (i = bset->n_ineq - 1; i >= 0; --i) {
 426                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + first, n) == -1)
 427                         continue;
 428                 isl_basic_set_drop_inequality(bset, i);
 429         }
 430
 431         return bset;
 432 }
 433
 434 /* Look for all equalities satisfied by the integer points in bset,
 435  * which is assumed to be bounded.
 436  *
 437  * The equalities are obtained by successively looking for
 438  * a point that is affinely independent of the points found so far.
 439  * In particular, for each equality satisfied by the points so far,
 440  * we check if there is any point on a hyperplane parallel to the
 441  * corresponding hyperplane shifted by at least one (in either direction).
 442  */
 443 static struct isl_basic_set *uset_affine_hull_bounded(struct isl_basic_set *bset)
 444 {
 445         struct isl_vec *sample = NULL;
 446         struct isl_basic_set *hull;
 447         struct isl_tab *tab = NULL;
 448         unsigned dim;
 449
 450         if (isl_basic_set_fast_is_empty(bset))
 451                 return bset;
 452
 453         dim = isl_basic_set_n_dim(bset);
 454
 455         if (bset->sample && bset->sample->size == 1 + dim) {
 456                 int contains = isl_basic_set_contains(bset, bset->sample);
 457                 if (contains < 0)
 458                         goto error;
 459                 if (contains) {
 460                         if (dim == 0)
 461                                 return bset;
 462                         sample = isl_vec_copy(bset->sample);
 463                 } else {
 464                         isl_vec_free(bset->sample);
 465                         bset->sample = NULL;
 466                 }
 467         }
 468
 469         tab = isl_tab_from_basic_set(bset);
 470         if (!tab)
 471                 goto error;
 472         if (tab->empty) {
 473                 isl_tab_free(tab);
 474                 isl_vec_free(sample);
 475                 return isl_basic_set_set_to_empty(bset);
 476         }
 477         if (isl_tab_track_bset(tab, isl_basic_set_copy(bset)) < 0)
 478                 goto error;
 479
 480         if (!sample) {
 481                 struct isl_tab_undo *snap;
 482                 snap = isl_tab_snap(tab);
 483                 sample = isl_tab_sample(tab);
 484                 if (isl_tab_rollback(tab, snap) < 0)
 485                         goto error;
 486                 isl_vec_free(tab->bmap->sample);
 487                 tab->bmap->sample = isl_vec_copy(sample);
 488         }
 489
 490         if (!sample)
 491                 goto error;
 492         if (sample->size == 0) {
 493                 isl_tab_free(tab);
 494                 isl_vec_free(sample);
 495                 return isl_basic_set_set_to_empty(bset);
 496         }
 497
 498         hull = isl_basic_set_from_vec(sample);
 499
 500         isl_basic_set_free(bset);
 501         hull = extend_affine_hull(tab, hull);
 502         isl_tab_free(tab);
 503
 504         return hull;
 505 error:
 506         isl_vec_free(sample);
 507         isl_tab_free(tab);
 508         isl_basic_set_free(bset);
 509         return NULL;
 510 }
 511
 512 /* Given an unbounded tableau and an integer point satisfying the tableau,
 513  * construct an intial affine hull containing the recession cone
 514  * shifted to the given point.
 515  *
 516  * The unbounded directions are taken from the last rows of the basis,
 517  * which is assumed to have been initialized appropriately.
 518  */
 519 static __isl_give isl_basic_set *initial_hull(struct isl_tab *tab,
 520         __isl_take isl_vec *vec)
 521 {
 522         int i;
 523         int k;
 524         struct isl_basic_set *bset = NULL;
 525         struct isl_ctx *ctx;
 526         unsigned dim;
 527
 528         if (!vec || !tab)
 529                 return NULL;
 530         ctx = vec->ctx;
 531         isl_assert(ctx, vec->size != 0, goto error);
 532
 533         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
 534         if (!bset)
 535                 goto error;
 536         dim = isl_basic_set_n_dim(bset) - tab->n_unbounded;
 537         for (i = 0; i < dim; ++i) {
 538                 k = isl_basic_set_alloc_equality(bset);
 539                 if (k < 0)
 540                         goto error;
 541                 isl_seq_cpy(bset->eq[k] + 1, tab->basis->row[1 + i] + 1,
 542                             vec->size - 1);
 543                 isl_seq_inner_product(bset->eq[k] + 1, vec->el +1,
 544                                       vec->size - 1, &bset->eq[k][0]);
 545                 isl_int_neg(bset->eq[k][0], bset->eq[k][0]);
 546         }
 547         bset->sample = vec;
 548         bset = isl_basic_set_gauss(bset, NULL);
 549
 550         return bset;
 551 error:
 552         isl_basic_set_free(bset);
 553         isl_vec_free(vec);
 554         return NULL;
 555 }
 556
 557 /* Given a tableau of a set and a tableau of the corresponding
 558  * recession cone, detect and add all equalities to the tableau.
 559  * If the tableau is bounded, then we can simply keep the
 560  * tableau in its state after the return from extend_affine_hull.
 561  * However, if the tableau is unbounded, then
 562  * isl_tab_set_initial_basis_with_cone will add some additional
 563  * constraints to the tableau that have to be removed again.
 564  * In this case, we therefore rollback to the state before
 565  * any constraints were added and then add the eqaulities back in.
 566  */
 567 struct isl_tab *isl_tab_detect_equalities(struct isl_tab *tab,
 568         struct isl_tab *tab_cone)
 569 {
 570         int j;
 571         struct isl_vec *sample;
 572         struct isl_basic_set *hull;
 573         struct isl_tab_undo *snap;
 574
 575         if (!tab || !tab_cone)
 576                 goto error;
 577
 578         snap = isl_tab_snap(tab);
 579
 580         isl_mat_free(tab->basis);
 581         tab->basis = NULL;
 582
 583         isl_assert(tab->mat->ctx, tab->bmap, goto error);
 584         isl_assert(tab->mat->ctx, tab->samples, goto error);
 585         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
 586         isl_assert(tab->mat->ctx, tab->n_sample > tab->n_outside, goto error);
 587
 588         if (isl_tab_set_initial_basis_with_cone(tab, tab_cone) < 0)
 589                 goto error;
 590
 591         sample = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
 592         if (!sample)
 593                 goto error;
 594
 595         isl_seq_cpy(sample->el, tab->samples->row[tab->n_outside], sample->size);
 596
 597         isl_vec_free(tab->bmap->sample);
 598         tab->bmap->sample = isl_vec_copy(sample);
 599
 600         if (tab->n_unbounded == 0)
 601                 hull = isl_basic_set_from_vec(isl_vec_copy(sample));
 602         else
 603                 hull = initial_hull(tab, isl_vec_copy(sample));
 604
 605         for (j = tab->n_outside + 1; j < tab->n_sample; ++j) {
 606                 isl_seq_cpy(sample->el, tab->samples->row[j], sample->size);
 607                 hull = affine_hull(hull,
 608                                 isl_basic_set_from_vec(isl_vec_copy(sample)));
 609         }
 610
 611         isl_vec_free(sample);
 612
 613         hull = extend_affine_hull(tab, hull);
 614         if (!hull)
 615                 goto error;
 616
 617         if (tab->n_unbounded == 0) {
 618                 isl_basic_set_free(hull);
 619                 return tab;
 620         }
 621
 622         if (isl_tab_rollback(tab, snap) < 0)
 623                 goto error;
 624
 625         if (hull->n_eq > tab->n_zero) {
 626                 for (j = 0; j < hull->n_eq; ++j) {
 627                         isl_seq_normalize(tab->mat->ctx, hull->eq[j], 1 + tab->n_var);
 628                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 629                 }
 630         }
 631
 632         isl_basic_set_free(hull);
 633
 634         return tab;
 635 error:
 636         isl_tab_free(tab);
 637         return NULL;
 638 }
 639
 640 /* Compute the affine hull of "bset", where "cone" is the recession cone
 641  * of "bset".
 642  *
 643  * We first compute a unimodular transformation that puts the unbounded
 644  * directions in the last dimensions.  In particular, we take a transformation
 645  * that maps all equalities to equalities (in HNF) on the first dimensions.
 646  * Let x be the original dimensions and y the transformed, with y_1 bounded
 647  * and y_2 unbounded.
 648  *
 649  *             [ y_1 ]                  [ y_1 ]   [ Q_1 ]
 650  *      x = U  [ y_2 ]                  [ y_2 ] = [ Q_2 ] x
 651  *
 652  * Let's call the input basic set S.  We compute S' = preimage(S, U)
 653  * and drop the final dimensions including any constraints involving them.
 654  * This results in set S''.
 655  * Then we compute the affine hull A'' of S''.
 656  * Let F y_1 >= g be the constraint system of A''.  In the transformed
 657  * space the y_2 are unbounded, so we can add them back without any constraints,
 658  * resulting in
 659  *
 660  *                      [ y_1 ]
 661  *              [ F 0 ] [ y_2 ] >= g
 662  * or
 663  *                      [ Q_1 ]
 664  *              [ F 0 ] [ Q_2 ] x >= g
 665  * or
 666  *              F Q_1 x >= g
 667  *
 668  * The affine hull in the original space is then obtained as
 669  * A = preimage(A'', Q_1).
 670  */
 671 static struct isl_basic_set *affine_hull_with_cone(struct isl_basic_set *bset,
 672         struct isl_basic_set *cone)
 673 {
 674         unsigned total;
 675         unsigned cone_dim;
 676         struct isl_basic_set *hull;
 677         struct isl_mat *M, *U, *Q;
 678
 679         if (!bset || !cone)
 680                 goto error;
 681
 682         total = isl_basic_set_total_dim(cone);
 683         cone_dim = total - cone->n_eq;
 684
 685         M = isl_mat_sub_alloc(bset->ctx, cone->eq, 0, cone->n_eq, 1, total);
 686         M = isl_mat_left_hermite(M, 0, &U, &Q);
 687         if (!M)
 688                 goto error;
 689         isl_mat_free(M);
 690
 691         U = isl_mat_lin_to_aff(U);
 692         bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
 693
 694         bset = drop_constraints_involving(bset, total - cone_dim, cone_dim);
 695         bset = isl_basic_set_drop_dims(bset, total - cone_dim, cone_dim);
 696
 697         Q = isl_mat_lin_to_aff(Q);
 698         Q = isl_mat_drop_rows(Q, 1 + total - cone_dim, cone_dim);
 699
 700         if (bset && bset->sample && bset->sample->size == 1 + total)
 701                 bset->sample = isl_mat_vec_product(isl_mat_copy(Q), bset->sample);
 702
 703         hull = uset_affine_hull_bounded(bset);
 704
 705         if (!hull)
 706                 isl_mat_free(U);
 707         else {
 708                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 709                 U = isl_mat_drop_cols(U, 1 + total - cone_dim, cone_dim);
 710                 if (sample && sample->size > 0)
 711                         sample = isl_mat_vec_product(U, sample);
 712                 else
 713                         isl_mat_free(U);
 714                 hull = isl_basic_set_preimage(hull, Q);
 715                 isl_vec_free(hull->sample);
 716                 hull->sample = sample;
 717         }
 718
 719         isl_basic_set_free(cone);
 720
 721         return hull;
 722 error:
 723         isl_basic_set_free(bset);
 724         isl_basic_set_free(cone);
 725         return NULL;
 726 }
 727
 728 /* Look for all equalities satisfied by the integer points in bset,
 729  * which is assumed not to have any explicit equalities.
 730  *
 731  * The equalities are obtained by successively looking for
 732  * a point that is affinely independent of the points found so far.
 733  * In particular, for each equality satisfied by the points so far,
 734  * we check if there is any point on a hyperplane parallel to the
 735  * corresponding hyperplane shifted by at least one (in either direction).
 736  *
 737  * Before looking for any outside points, we first compute the recession
 738  * cone.  The directions of this recession cone will always be part
 739  * of the affine hull, so there is no need for looking for any points
 740  * in these directions.
 741  * In particular, if the recession cone is full-dimensional, then
 742  * the affine hull is simply the whole universe.
 743  */
 744 static struct isl_basic_set *uset_affine_hull(struct isl_basic_set *bset)
 745 {
 746         struct isl_basic_set *cone;
 747
 748         if (isl_basic_set_fast_is_empty(bset))
 749                 return bset;
 750
 751         cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
 752         if (!cone)
 753                 goto error;
 754         if (cone->n_eq == 0) {
 755                 struct isl_basic_set *hull;
 756                 isl_basic_set_free(cone);
 757                 hull = isl_basic_set_universe_like(bset);
 758                 isl_basic_set_free(bset);
 759                 return hull;
 760         }
 761
 762         if (cone->n_eq < isl_basic_set_total_dim(cone))
 763                 return affine_hull_with_cone(bset, cone);
 764
 765         isl_basic_set_free(cone);
 766         return uset_affine_hull_bounded(bset);
 767 error:
 768         isl_basic_set_free(bset);
 769         return NULL;
 770 }
 771
 772 /* Look for all equalities satisfied by the integer points in bmap
 773  * that are independent of the equalities already explicitly available
 774  * in bmap.
 775  *
 776  * We first remove all equalities already explicitly available,
 777  * then look for additional equalities in the reduced space
 778  * and then transform the result to the original space.
 779  * The original equalities are _not_ added to this set.  This is
 780  * the responsibility of the calling function.
 781  * The resulting basic set has all meaning about the dimensions removed.
 782  * In particular, dimensions that correspond to existential variables
 783  * in bmap and that are found to be fixed are not removed.
 784  */
 785 static struct isl_basic_set *equalities_in_underlying_set(
 786                                                 struct isl_basic_map *bmap)
 787 {
 788         struct isl_mat *T1 = NULL;
 789         struct isl_mat *T2 = NULL;
 790         struct isl_basic_set *bset = NULL;
 791         struct isl_basic_set *hull = NULL;
 792
 793         bset = isl_basic_map_underlying_set(bmap);
 794         if (!bset)
 795                 return NULL;
 796         if (bset->n_eq)
 797                 bset = isl_basic_set_remove_equalities(bset, &T1, &T2);
 798         if (!bset)
 799                 goto error;
 800
 801         hull = uset_affine_hull(bset);
 802         if (!T2)
 803                 return hull;
 804
 805         if (!hull)
 806                 isl_mat_free(T1);
 807         else {
 808                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 809                 if (sample && sample->size > 0)
 810                         sample = isl_mat_vec_product(T1, sample);
 811                 else
 812                         isl_mat_free(T1);
 813                 hull = isl_basic_set_preimage(hull, T2);
 814                 isl_vec_free(hull->sample);
 815                 hull->sample = sample;
 816         }
 817
 818         return hull;
 819 error:
 820         isl_mat_free(T2);
 821         isl_basic_set_free(bset);
 822         isl_basic_set_free(hull);
 823         return NULL;
 824 }
 825
 826 /* Detect and make explicit all equalities satisfied by the (integer)
 827  * points in bmap.
 828  */
 829 struct isl_basic_map *isl_basic_map_detect_equalities(
 830                                                 struct isl_basic_map *bmap)
 831 {
 832         int i, j;
 833         struct isl_basic_set *hull = NULL;
 834
 835         if (!bmap)
 836                 return NULL;
 837         if (bmap->n_ineq == 0)
 838                 return bmap;
 839         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
 840                 return bmap;
 841         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_ALL_EQUALITIES))
 842                 return bmap;
 843         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 844                 return isl_basic_map_implicit_equalities(bmap);
 845
 846         hull = equalities_in_underlying_set(isl_basic_map_copy(bmap));
 847         if (!hull)
 848                 goto error;
 849         if (ISL_F_ISSET(hull, ISL_BASIC_SET_EMPTY)) {
 850                 isl_basic_set_free(hull);
 851                 return isl_basic_map_set_to_empty(bmap);
 852         }
 853         bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim), 0,
 854                                         hull->n_eq, 0);
 855         for (i = 0; i < hull->n_eq; ++i) {
 856                 j = isl_basic_map_alloc_equality(bmap);
 857                 if (j < 0)
 858                         goto error;
 859                 isl_seq_cpy(bmap->eq[j], hull->eq[i],
 860                                 1 + isl_basic_set_total_dim(hull));
 861         }
 862         isl_vec_free(bmap->sample);
 863         bmap->sample = isl_vec_copy(hull->sample);
 864         isl_basic_set_free(hull);
 865         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT | ISL_BASIC_MAP_ALL_EQUALITIES);
 866         bmap = isl_basic_map_simplify(bmap);
 867         return isl_basic_map_finalize(bmap);
 868 error:
 869         isl_basic_set_free(hull);
 870         isl_basic_map_free(bmap);
 871         return NULL;
 872 }
 873
 874 __isl_give isl_basic_set *isl_basic_set_detect_equalities(
 875                                                 __isl_take isl_basic_set *bset)
 876 {
 877         return (isl_basic_set *)
 878                 isl_basic_map_detect_equalities((isl_basic_map *)bset);
 879 }
 880
 881 struct isl_map *isl_map_detect_equalities(struct isl_map *map)
 882 {
 883         struct isl_basic_map *bmap;
 884         int i;
 885
 886         if (!map)
 887                 return NULL;
 888
 889         for (i = 0; i < map->n; ++i) {
 890                 bmap = isl_basic_map_copy(map->p[i]);
 891                 bmap = isl_basic_map_detect_equalities(bmap);
 892                 if (!bmap)
 893                         goto error;
 894                 isl_basic_map_free(map->p[i]);
 895                 map->p[i] = bmap;
 896         }
 897
 898         return map;
 899 error:
 900         isl_map_free(map);
 901         return NULL;
 902 }
 903
 904 __isl_give isl_set *isl_set_detect_equalities(__isl_take isl_set *set)
 905 {
 906         return (isl_set *)isl_map_detect_equalities((isl_map *)set);
 907 }
 908
 909 /* After computing the rational affine hull (by detecting the implicit
 910  * equalities), we compute the additional equalities satisfied by
 911  * the integer points (if any) and add the original equalities back in.
 912  */
 913 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_basic_map *bmap)
 914 {
 915         bmap = isl_basic_map_detect_equalities(bmap);
 916         bmap = isl_basic_map_cow(bmap);
 917         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 918         return bmap;
 919 }
 920
 921 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_basic_set *bset)
 922 {
 923         return (struct isl_basic_set *)
 924                 isl_basic_map_affine_hull((struct isl_basic_map *)bset);
 925 }
 926
 927 struct isl_basic_map *isl_map_affine_hull(struct isl_map *map)
 928 {
 929         int i;
 930         struct isl_basic_map *model = NULL;
 931         struct isl_basic_map *hull = NULL;
 932         struct isl_set *set;
 933
 934         map = isl_map_detect_equalities(map);
 935         map = isl_map_align_divs(map);
 936
 937         if (!map)
 938                 return NULL;
 939
 940         if (map->n == 0) {
 941                 hull = isl_basic_map_empty_like_map(map);
 942                 isl_map_free(map);
 943                 return hull;
 944         }
 945
 946         model = isl_basic_map_copy(map->p[0]);
 947         set = isl_map_underlying_set(map);
 948         set = isl_set_cow(set);
 949         if (!set)
 950                 goto error;
 951
 952         for (i = 0; i < set->n; ++i) {
 953                 set->p[i] = isl_basic_set_cow(set->p[i]);
 954                 set->p[i] = isl_basic_set_affine_hull(set->p[i]);
 955                 set->p[i] = isl_basic_set_gauss(set->p[i], NULL);
 956                 if (!set->p[i])
 957                         goto error;
 958         }
 959         set = isl_set_remove_empty_parts(set);
 960         if (set->n == 0) {
 961                 hull = isl_basic_map_empty_like(model);
 962                 isl_basic_map_free(model);
 963         } else {
 964                 struct isl_basic_set *bset;
 965                 while (set->n > 1) {
 966                         set->p[0] = affine_hull(set->p[0], set->p[--set->n]);
 967                         if (!set->p[0])
 968                                 goto error;
 969                 }
 970                 bset = isl_basic_set_copy(set->p[0]);
 971                 hull = isl_basic_map_overlying_set(bset, model);
 972         }
 973         isl_set_free(set);
 974         hull = isl_basic_map_simplify(hull);
 975         return isl_basic_map_finalize(hull);
 976 error:
 977         isl_basic_map_free(model);
 978         isl_set_free(set);
 979         return NULL;
 980 }
 981
 982 struct isl_basic_set *isl_set_affine_hull(struct isl_set *set)
 983 {
 984         return (struct isl_basic_set *)
 985                 isl_map_affine_hull((struct isl_map *)set);
 986 }