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         int k;
 255
 256         if (!tab)
 257                 return NULL;
 258         ctx = tab->mat->ctx;
 259
 260         dim = tab->n_var;
 261         sample = isl_vec_alloc(ctx, 1 + dim);
 262         if (!sample)
 263                 return NULL;
 264         isl_int_set_si(sample->el[0], 1);
 265         isl_seq_combine(sample->el + 1,
 266                 ctx->one, tab->bmap->sample->el + 1,
 267                 up ? ctx->one : ctx->negone, eq + 1, dim);
 268         if (isl_basic_map_contains(tab->bmap, sample))
 269                 return sample;
 270         isl_vec_free(sample);
 271         sample = NULL;
 272
 273         snap = isl_tab_snap(tab);
 274
 275         if (!up)
 276                 isl_seq_neg(eq, eq, 1 + dim);
 277         isl_int_sub_ui(eq[0], eq[0], 1);
 278
 279         if (isl_tab_extend_cons(tab, 1) < 0)
 280                 goto error;
 281         if (isl_tab_add_ineq(tab, eq) < 0)
 282                 goto error;
 283
 284         sample = isl_tab_sample(tab);
 285
 286         isl_int_add_ui(eq[0], eq[0], 1);
 287         if (!up)
 288                 isl_seq_neg(eq, eq, 1 + dim);
 289
 290         if (isl_tab_rollback(tab, snap) < 0)
 291                 goto error;
 292
 293         return sample;
 294 error:
 295         isl_vec_free(sample);
 296         return NULL;
 297 }
 298
 299 struct isl_basic_set *isl_basic_set_recession_cone(struct isl_basic_set *bset)
 300 {
 301         int i;
 302
 303         bset = isl_basic_set_cow(bset);
 304         if (!bset)
 305                 return NULL;
 306         isl_assert(bset->ctx, bset->n_div == 0, goto error);
 307
 308         for (i = 0; i < bset->n_eq; ++i)
 309                 isl_int_set_si(bset->eq[i][0], 0);
 310
 311         for (i = 0; i < bset->n_ineq; ++i)
 312                 isl_int_set_si(bset->ineq[i][0], 0);
 313
 314         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
 315         return isl_basic_set_implicit_equalities(bset);
 316 error:
 317         isl_basic_set_free(bset);
 318         return NULL;
 319 }
 320
 321 __isl_give isl_set *isl_set_recession_cone(__isl_take isl_set *set)
 322 {
 323         int i;
 324
 325         if (!set)
 326                 return NULL;
 327         if (set->n == 0)
 328                 return set;
 329
 330         set = isl_set_remove_divs(set);
 331         set = isl_set_cow(set);
 332         if (!set)
 333                 return NULL;
 334
 335         for (i = 0; i < set->n; ++i) {
 336                 set->p[i] = isl_basic_set_recession_cone(set->p[i]);
 337                 if (!set->p[i])
 338                         goto error;
 339         }
 340
 341         return set;
 342 error:
 343         isl_set_free(set);
 344         return NULL;
 345 }
 346
 347 /* Extend an initial (under-)approximation of the affine hull of basic
 348  * set represented by the tableau "tab"
 349  * by looking for points that do not satisfy one of the equalities
 350  * in the current approximation and adding them to that approximation
 351  * until no such points can be found any more.
 352  *
 353  * The caller of this function ensures that "tab" is bounded or
 354  * that tab->basis and tab->n_unbounded have been set appropriately.
 355  */
 356 static struct isl_basic_set *extend_affine_hull(struct isl_tab *tab,
 357         struct isl_basic_set *hull)
 358 {
 359         int i, j, k;
 360         unsigned dim;
 361
 362         if (!tab || !hull)
 363                 goto error;
 364
 365         dim = tab->n_var;
 366
 367         if (isl_tab_extend_cons(tab, 2 * dim + 1) < 0)
 368                 goto error;
 369
 370         for (i = 0; i < dim; ++i) {
 371                 struct isl_vec *sample;
 372                 struct isl_basic_set *point;
 373                 for (j = 0; j < hull->n_eq; ++j) {
 374                         sample = outside_point(tab, hull->eq[j], 1);
 375                         if (!sample)
 376                                 goto error;
 377                         if (sample->size > 0)
 378                                 break;
 379                         isl_vec_free(sample);
 380                         sample = outside_point(tab, hull->eq[j], 0);
 381                         if (!sample)
 382                                 goto error;
 383                         if (sample->size > 0)
 384                                 break;
 385                         isl_vec_free(sample);
 386
 387                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 388                         if (!tab)
 389                                 goto error;
 390                 }
 391                 if (j == hull->n_eq)
 392                         break;
 393                 if (tab->samples)
 394                         tab = isl_tab_add_sample(tab, isl_vec_copy(sample));
 395                 if (!tab)
 396                         goto error;
 397                 point = isl_basic_set_from_vec(sample);
 398                 hull = affine_hull(hull, point);
 399         }
 400
 401         return hull;
 402 error:
 403         isl_basic_set_free(hull);
 404         return NULL;
 405 }
 406
 407 /* Drop all constraints in bset that involve any of the dimensions
 408  * first to first+n-1.
 409  */
 410 static struct isl_basic_set *drop_constraints_involving
 411         (struct isl_basic_set *bset, unsigned first, unsigned n)
 412 {
 413         int i;
 414
 415         if (!bset)
 416                 return NULL;
 417
 418         bset = isl_basic_set_cow(bset);
 419
 420         for (i = bset->n_eq - 1; i >= 0; --i) {
 421                 if (isl_seq_first_non_zero(bset->eq[i] + 1 + first, n) == -1)
 422                         continue;
 423                 isl_basic_set_drop_equality(bset, i);
 424         }
 425
 426         for (i = bset->n_ineq - 1; i >= 0; --i) {
 427                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + first, n) == -1)
 428                         continue;
 429                 isl_basic_set_drop_inequality(bset, i);
 430         }
 431
 432         return bset;
 433 }
 434
 435 /* Look for all equalities satisfied by the integer points in bset,
 436  * which is assumed to be bounded.
 437  *
 438  * The equalities are obtained by successively looking for
 439  * a point that is affinely independent of the points found so far.
 440  * In particular, for each equality satisfied by the points so far,
 441  * we check if there is any point on a hyperplane parallel to the
 442  * corresponding hyperplane shifted by at least one (in either direction).
 443  */
 444 static struct isl_basic_set *uset_affine_hull_bounded(struct isl_basic_set *bset)
 445 {
 446         struct isl_vec *sample = NULL;
 447         struct isl_basic_set *hull;
 448         struct isl_tab *tab = NULL;
 449         unsigned dim;
 450
 451         if (isl_basic_set_fast_is_empty(bset))
 452                 return bset;
 453
 454         dim = isl_basic_set_n_dim(bset);
 455
 456         if (bset->sample && bset->sample->size == 1 + dim) {
 457                 int contains = isl_basic_set_contains(bset, bset->sample);
 458                 if (contains < 0)
 459                         goto error;
 460                 if (contains) {
 461                         if (dim == 0)
 462                                 return bset;
 463                         sample = isl_vec_copy(bset->sample);
 464                 } else {
 465                         isl_vec_free(bset->sample);
 466                         bset->sample = NULL;
 467                 }
 468         }
 469
 470         tab = isl_tab_from_basic_set(bset);
 471         if (!tab)
 472                 goto error;
 473         if (tab->empty) {
 474                 isl_tab_free(tab);
 475                 isl_vec_free(sample);
 476                 return isl_basic_set_set_to_empty(bset);
 477         }
 478         if (isl_tab_track_bset(tab, isl_basic_set_copy(bset)) < 0)
 479                 goto error;
 480
 481         if (!sample) {
 482                 struct isl_tab_undo *snap;
 483                 snap = isl_tab_snap(tab);
 484                 sample = isl_tab_sample(tab);
 485                 if (isl_tab_rollback(tab, snap) < 0)
 486                         goto error;
 487                 isl_vec_free(tab->bmap->sample);
 488                 tab->bmap->sample = isl_vec_copy(sample);
 489         }
 490
 491         if (!sample)
 492                 goto error;
 493         if (sample->size == 0) {
 494                 isl_tab_free(tab);
 495                 isl_vec_free(sample);
 496                 return isl_basic_set_set_to_empty(bset);
 497         }
 498
 499         hull = isl_basic_set_from_vec(sample);
 500
 501         isl_basic_set_free(bset);
 502         hull = extend_affine_hull(tab, hull);
 503         isl_tab_free(tab);
 504
 505         return hull;
 506 error:
 507         isl_vec_free(sample);
 508         isl_tab_free(tab);
 509         isl_basic_set_free(bset);
 510         return NULL;
 511 }
 512
 513 /* Given an unbounded tableau and an integer point satisfying the tableau,
 514  * construct an intial affine hull containing the recession cone
 515  * shifted to the given point.
 516  *
 517  * The unbounded directions are taken from the last rows of the basis,
 518  * which is assumed to have been initialized appropriately.
 519  */
 520 static __isl_give isl_basic_set *initial_hull(struct isl_tab *tab,
 521         __isl_take isl_vec *vec)
 522 {
 523         int i;
 524         int k;
 525         struct isl_basic_set *bset = NULL;
 526         struct isl_ctx *ctx;
 527         unsigned dim;
 528
 529         if (!vec || !tab)
 530                 return NULL;
 531         ctx = vec->ctx;
 532         isl_assert(ctx, vec->size != 0, goto error);
 533
 534         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
 535         if (!bset)
 536                 goto error;
 537         dim = isl_basic_set_n_dim(bset) - tab->n_unbounded;
 538         for (i = 0; i < dim; ++i) {
 539                 k = isl_basic_set_alloc_equality(bset);
 540                 if (k < 0)
 541                         goto error;
 542                 isl_seq_cpy(bset->eq[k] + 1, tab->basis->row[1 + i] + 1,
 543                             vec->size - 1);
 544                 isl_seq_inner_product(bset->eq[k] + 1, vec->el +1,
 545                                       vec->size - 1, &bset->eq[k][0]);
 546                 isl_int_neg(bset->eq[k][0], bset->eq[k][0]);
 547         }
 548         bset->sample = vec;
 549         bset = isl_basic_set_gauss(bset, NULL);
 550
 551         return bset;
 552 error:
 553         isl_basic_set_free(bset);
 554         isl_vec_free(vec);
 555         return NULL;
 556 }
 557
 558 /* Given a tableau of a set and a tableau of the corresponding
 559  * recession cone, detect and add all equalities to the tableau.
 560  * If the tableau is bounded, then we can simply keep the
 561  * tableau in its state after the return from extend_affine_hull.
 562  * However, if the tableau is unbounded, then
 563  * isl_tab_set_initial_basis_with_cone will add some additional
 564  * constraints to the tableau that have to be removed again.
 565  * In this case, we therefore rollback to the state before
 566  * any constraints were added and then add the eqaulities back in.
 567  */
 568 struct isl_tab *isl_tab_detect_equalities(struct isl_tab *tab,
 569         struct isl_tab *tab_cone)
 570 {
 571         int j;
 572         struct isl_vec *sample;
 573         struct isl_basic_set *hull;
 574         struct isl_tab_undo *snap;
 575
 576         if (!tab || !tab_cone)
 577                 goto error;
 578
 579         snap = isl_tab_snap(tab);
 580
 581         isl_mat_free(tab->basis);
 582         tab->basis = NULL;
 583
 584         isl_assert(tab->mat->ctx, tab->bmap, goto error);
 585         isl_assert(tab->mat->ctx, tab->samples, goto error);
 586         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
 587         isl_assert(tab->mat->ctx, tab->n_sample > tab->n_outside, goto error);
 588
 589         if (isl_tab_set_initial_basis_with_cone(tab, tab_cone) < 0)
 590                 goto error;
 591
 592         sample = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
 593         if (!sample)
 594                 goto error;
 595
 596         isl_seq_cpy(sample->el, tab->samples->row[tab->n_outside], sample->size);
 597
 598         isl_vec_free(tab->bmap->sample);
 599         tab->bmap->sample = isl_vec_copy(sample);
 600
 601         if (tab->n_unbounded == 0)
 602                 hull = isl_basic_set_from_vec(isl_vec_copy(sample));
 603         else
 604                 hull = initial_hull(tab, isl_vec_copy(sample));
 605
 606         for (j = tab->n_outside + 1; j < tab->n_sample; ++j) {
 607                 isl_seq_cpy(sample->el, tab->samples->row[j], sample->size);
 608                 hull = affine_hull(hull,
 609                                 isl_basic_set_from_vec(isl_vec_copy(sample)));
 610         }
 611
 612         isl_vec_free(sample);
 613
 614         hull = extend_affine_hull(tab, hull);
 615         if (!hull)
 616                 goto error;
 617
 618         if (tab->n_unbounded == 0) {
 619                 isl_basic_set_free(hull);
 620                 return tab;
 621         }
 622
 623         if (isl_tab_rollback(tab, snap) < 0)
 624                 goto error;
 625
 626         if (hull->n_eq > tab->n_zero) {
 627                 for (j = 0; j < hull->n_eq; ++j) {
 628                         isl_seq_normalize(tab->mat->ctx, hull->eq[j], 1 + tab->n_var);
 629                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 630                 }
 631         }
 632
 633         isl_basic_set_free(hull);
 634
 635         return tab;
 636 error:
 637         isl_tab_free(tab);
 638         return NULL;
 639 }
 640
 641 /* Compute the affine hull of "bset", where "cone" is the recession cone
 642  * of "bset".
 643  *
 644  * We first compute a unimodular transformation that puts the unbounded
 645  * directions in the last dimensions.  In particular, we take a transformation
 646  * that maps all equalities to equalities (in HNF) on the first dimensions.
 647  * Let x be the original dimensions and y the transformed, with y_1 bounded
 648  * and y_2 unbounded.
 649  *
 650  *             [ y_1 ]                  [ y_1 ]   [ Q_1 ]
 651  *      x = U  [ y_2 ]                  [ y_2 ] = [ Q_2 ] x
 652  *
 653  * Let's call the input basic set S.  We compute S' = preimage(S, U)
 654  * and drop the final dimensions including any constraints involving them.
 655  * This results in set S''.
 656  * Then we compute the affine hull A'' of S''.
 657  * Let F y_1 >= g be the constraint system of A''.  In the transformed
 658  * space the y_2 are unbounded, so we can add them back without any constraints,
 659  * resulting in
 660  *
 661  *                      [ y_1 ]
 662  *              [ F 0 ] [ y_2 ] >= g
 663  * or
 664  *                      [ Q_1 ]
 665  *              [ F 0 ] [ Q_2 ] x >= g
 666  * or
 667  *              F Q_1 x >= g
 668  *
 669  * The affine hull in the original space is then obtained as
 670  * A = preimage(A'', Q_1).
 671  */
 672 static struct isl_basic_set *affine_hull_with_cone(struct isl_basic_set *bset,
 673         struct isl_basic_set *cone)
 674 {
 675         unsigned total;
 676         unsigned cone_dim;
 677         struct isl_basic_set *hull;
 678         struct isl_mat *M, *U, *Q;
 679
 680         if (!bset || !cone)
 681                 goto error;
 682
 683         total = isl_basic_set_total_dim(cone);
 684         cone_dim = total - cone->n_eq;
 685
 686         M = isl_mat_sub_alloc(bset->ctx, cone->eq, 0, cone->n_eq, 1, total);
 687         M = isl_mat_left_hermite(M, 0, &U, &Q);
 688         if (!M)
 689                 goto error;
 690         isl_mat_free(M);
 691
 692         U = isl_mat_lin_to_aff(U);
 693         bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
 694
 695         bset = drop_constraints_involving(bset, total - cone_dim, cone_dim);
 696         bset = isl_basic_set_drop_dims(bset, total - cone_dim, cone_dim);
 697
 698         Q = isl_mat_lin_to_aff(Q);
 699         Q = isl_mat_drop_rows(Q, 1 + total - cone_dim, cone_dim);
 700
 701         if (bset && bset->sample && bset->sample->size == 1 + total)
 702                 bset->sample = isl_mat_vec_product(isl_mat_copy(Q), bset->sample);
 703
 704         hull = uset_affine_hull_bounded(bset);
 705
 706         if (!hull)
 707                 isl_mat_free(U);
 708         else {
 709                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 710                 U = isl_mat_drop_cols(U, 1 + total - cone_dim, cone_dim);
 711                 if (sample && sample->size > 0)
 712                         sample = isl_mat_vec_product(U, sample);
 713                 else
 714                         isl_mat_free(U);
 715                 hull = isl_basic_set_preimage(hull, Q);
 716                 isl_vec_free(hull->sample);
 717                 hull->sample = sample;
 718         }
 719
 720         isl_basic_set_free(cone);
 721
 722         return hull;
 723 error:
 724         isl_basic_set_free(bset);
 725         isl_basic_set_free(cone);
 726         return NULL;
 727 }
 728
 729 /* Look for all equalities satisfied by the integer points in bset,
 730  * which is assumed not to have any explicit equalities.
 731  *
 732  * The equalities are obtained by successively looking for
 733  * a point that is affinely independent of the points found so far.
 734  * In particular, for each equality satisfied by the points so far,
 735  * we check if there is any point on a hyperplane parallel to the
 736  * corresponding hyperplane shifted by at least one (in either direction).
 737  *
 738  * Before looking for any outside points, we first compute the recession
 739  * cone.  The directions of this recession cone will always be part
 740  * of the affine hull, so there is no need for looking for any points
 741  * in these directions.
 742  * In particular, if the recession cone is full-dimensional, then
 743  * the affine hull is simply the whole universe.
 744  */
 745 static struct isl_basic_set *uset_affine_hull(struct isl_basic_set *bset)
 746 {
 747         struct isl_basic_set *cone;
 748
 749         if (isl_basic_set_fast_is_empty(bset))
 750                 return bset;
 751
 752         cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
 753         if (!cone)
 754                 goto error;
 755         if (cone->n_eq == 0) {
 756                 struct isl_basic_set *hull;
 757                 isl_basic_set_free(cone);
 758                 hull = isl_basic_set_universe_like(bset);
 759                 isl_basic_set_free(bset);
 760                 return hull;
 761         }
 762
 763         if (cone->n_eq < isl_basic_set_total_dim(cone))
 764                 return affine_hull_with_cone(bset, cone);
 765
 766         isl_basic_set_free(cone);
 767         return uset_affine_hull_bounded(bset);
 768 error:
 769         isl_basic_set_free(bset);
 770         return NULL;
 771 }
 772
 773 /* Look for all equalities satisfied by the integer points in bmap
 774  * that are independent of the equalities already explicitly available
 775  * in bmap.
 776  *
 777  * We first remove all equalities already explicitly available,
 778  * then look for additional equalities in the reduced space
 779  * and then transform the result to the original space.
 780  * The original equalities are _not_ added to this set.  This is
 781  * the responsibility of the calling function.
 782  * The resulting basic set has all meaning about the dimensions removed.
 783  * In particular, dimensions that correspond to existential variables
 784  * in bmap and that are found to be fixed are not removed.
 785  */
 786 static struct isl_basic_set *equalities_in_underlying_set(
 787                                                 struct isl_basic_map *bmap)
 788 {
 789         struct isl_mat *T1 = NULL;
 790         struct isl_mat *T2 = NULL;
 791         struct isl_basic_set *bset = NULL;
 792         struct isl_basic_set *hull = NULL;
 793
 794         bset = isl_basic_map_underlying_set(bmap);
 795         if (!bset)
 796                 return NULL;
 797         if (bset->n_eq)
 798                 bset = isl_basic_set_remove_equalities(bset, &T1, &T2);
 799         if (!bset)
 800                 goto error;
 801
 802         hull = uset_affine_hull(bset);
 803         if (!T2)
 804                 return hull;
 805
 806         if (!hull)
 807                 isl_mat_free(T1);
 808         else {
 809                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 810                 if (sample && sample->size > 0)
 811                         sample = isl_mat_vec_product(T1, sample);
 812                 else
 813                         isl_mat_free(T1);
 814                 hull = isl_basic_set_preimage(hull, T2);
 815                 isl_vec_free(hull->sample);
 816                 hull->sample = sample;
 817         }
 818
 819         return hull;
 820 error:
 821         isl_mat_free(T2);
 822         isl_basic_set_free(bset);
 823         isl_basic_set_free(hull);
 824         return NULL;
 825 }
 826
 827 /* Detect and make explicit all equalities satisfied by the (integer)
 828  * points in bmap.
 829  */
 830 struct isl_basic_map *isl_basic_map_detect_equalities(
 831                                                 struct isl_basic_map *bmap)
 832 {
 833         int i, j;
 834         struct isl_basic_set *hull = NULL;
 835
 836         if (!bmap)
 837                 return NULL;
 838         if (bmap->n_ineq == 0)
 839                 return bmap;
 840         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
 841                 return bmap;
 842         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_ALL_EQUALITIES))
 843                 return bmap;
 844         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 845                 return isl_basic_map_implicit_equalities(bmap);
 846
 847         hull = equalities_in_underlying_set(isl_basic_map_copy(bmap));
 848         if (!hull)
 849                 goto error;
 850         if (ISL_F_ISSET(hull, ISL_BASIC_SET_EMPTY)) {
 851                 isl_basic_set_free(hull);
 852                 return isl_basic_map_set_to_empty(bmap);
 853         }
 854         bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim), 0,
 855                                         hull->n_eq, 0);
 856         for (i = 0; i < hull->n_eq; ++i) {
 857                 j = isl_basic_map_alloc_equality(bmap);
 858                 if (j < 0)
 859                         goto error;
 860                 isl_seq_cpy(bmap->eq[j], hull->eq[i],
 861                                 1 + isl_basic_set_total_dim(hull));
 862         }
 863         isl_vec_free(bmap->sample);
 864         bmap->sample = isl_vec_copy(hull->sample);
 865         isl_basic_set_free(hull);
 866         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT | ISL_BASIC_MAP_ALL_EQUALITIES);
 867         bmap = isl_basic_map_simplify(bmap);
 868         return isl_basic_map_finalize(bmap);
 869 error:
 870         isl_basic_set_free(hull);
 871         isl_basic_map_free(bmap);
 872         return NULL;
 873 }
 874
 875 __isl_give isl_basic_set *isl_basic_set_detect_equalities(
 876                                                 __isl_take isl_basic_set *bset)
 877 {
 878         return (isl_basic_set *)
 879                 isl_basic_map_detect_equalities((isl_basic_map *)bset);
 880 }
 881
 882 struct isl_map *isl_map_detect_equalities(struct isl_map *map)
 883 {
 884         struct isl_basic_map *bmap;
 885         int i;
 886
 887         if (!map)
 888                 return NULL;
 889
 890         for (i = 0; i < map->n; ++i) {
 891                 bmap = isl_basic_map_copy(map->p[i]);
 892                 bmap = isl_basic_map_detect_equalities(bmap);
 893                 if (!bmap)
 894                         goto error;
 895                 isl_basic_map_free(map->p[i]);
 896                 map->p[i] = bmap;
 897         }
 898
 899         return map;
 900 error:
 901         isl_map_free(map);
 902         return NULL;
 903 }
 904
 905 __isl_give isl_set *isl_set_detect_equalities(__isl_take isl_set *set)
 906 {
 907         return (isl_set *)isl_map_detect_equalities((isl_map *)set);
 908 }
 909
 910 /* After computing the rational affine hull (by detecting the implicit
 911  * equalities), we compute the additional equalities satisfied by
 912  * the integer points (if any) and add the original equalities back in.
 913  */
 914 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_basic_map *bmap)
 915 {
 916         bmap = isl_basic_map_detect_equalities(bmap);
 917         bmap = isl_basic_map_cow(bmap);
 918         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 919         return bmap;
 920 }
 921
 922 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_basic_set *bset)
 923 {
 924         return (struct isl_basic_set *)
 925                 isl_basic_map_affine_hull((struct isl_basic_map *)bset);
 926 }
 927
 928 struct isl_basic_map *isl_map_affine_hull(struct isl_map *map)
 929 {
 930         int i;
 931         struct isl_basic_map *model = NULL;
 932         struct isl_basic_map *hull = NULL;
 933         struct isl_set *set;
 934
 935         if (!map)
 936                 return NULL;
 937
 938         if (map->n == 0) {
 939                 hull = isl_basic_map_empty_like_map(map);
 940                 isl_map_free(map);
 941                 return hull;
 942         }
 943
 944         map = isl_map_detect_equalities(map);
 945         map = isl_map_align_divs(map);
 946         if (!map)
 947                 return NULL;
 948         model = isl_basic_map_copy(map->p[0]);
 949         set = isl_map_underlying_set(map);
 950         set = isl_set_cow(set);
 951         if (!set)
 952                 goto error;
 953
 954         for (i = 0; i < set->n; ++i) {
 955                 set->p[i] = isl_basic_set_cow(set->p[i]);
 956                 set->p[i] = isl_basic_set_affine_hull(set->p[i]);
 957                 set->p[i] = isl_basic_set_gauss(set->p[i], NULL);
 958                 if (!set->p[i])
 959                         goto error;
 960         }
 961         set = isl_set_remove_empty_parts(set);
 962         if (set->n == 0) {
 963                 hull = isl_basic_map_empty_like(model);
 964                 isl_basic_map_free(model);
 965         } else {
 966                 struct isl_basic_set *bset;
 967                 while (set->n > 1) {
 968                         set->p[0] = affine_hull(set->p[0], set->p[--set->n]);
 969                         if (!set->p[0])
 970                                 goto error;
 971                 }
 972                 bset = isl_basic_set_copy(set->p[0]);
 973                 hull = isl_basic_map_overlying_set(bset, model);
 974         }
 975         isl_set_free(set);
 976         hull = isl_basic_map_simplify(hull);
 977         return isl_basic_map_finalize(hull);
 978 error:
 979         isl_basic_map_free(model);
 980         isl_set_free(set);
 981         return NULL;
 982 }
 983
 984 struct isl_basic_set *isl_set_affine_hull(struct isl_set *set)
 985 {
 986         return (struct isl_basic_set *)
 987                 isl_map_affine_hull((struct isl_map *)set);
 988 }