isl_sample.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the MIT 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_private.h>
  11 #include <isl_map_private.h>
  12 #include "isl_sample.h"
  13 #include <isl/vec.h>
  14 #include <isl/mat.h>
  15 #include <isl_seq.h>
  16 #include "isl_equalities.h"
  17 #include "isl_tab.h"
  18 #include "isl_basis_reduction.h"
  19 #include <isl_factorization.h>
  20 #include <isl_point_private.h>
  21 #include <isl_options_private.h>
  22 #include <isl_vec_private.h>
  23
  24 #include <bset_from_bmap.c>
  25 #include <set_to_map.c>
  26
  27 static __isl_give isl_vec *empty_sample(__isl_take isl_basic_set *bset)
  28 {
  29         struct isl_vec *vec;
  30
  31         vec = isl_vec_alloc(bset->ctx, 0);
  32         isl_basic_set_free(bset);
  33         return vec;
  34 }
  35
  36 /* Construct a zero sample of the same dimension as bset.
  37  * As a special case, if bset is zero-dimensional, this
  38  * function creates a zero-dimensional sample point.
  39  */
  40 static __isl_give isl_vec *zero_sample(__isl_take isl_basic_set *bset)
  41 {
  42         unsigned dim;
  43         struct isl_vec *sample;
  44
  45         dim = isl_basic_set_total_dim(bset);
  46         sample = isl_vec_alloc(bset->ctx, 1 + dim);
  47         if (sample) {
  48                 isl_int_set_si(sample->el[0], 1);
  49                 isl_seq_clr(sample->el + 1, dim);
  50         }
  51         isl_basic_set_free(bset);
  52         return sample;
  53 }
  54
  55 static __isl_give isl_vec *interval_sample(__isl_take isl_basic_set *bset)
  56 {
  57         int i;
  58         isl_int t;
  59         struct isl_vec *sample;
  60
  61         bset = isl_basic_set_simplify(bset);
  62         if (!bset)
  63                 return NULL;
  64         if (isl_basic_set_plain_is_empty(bset))
  65                 return empty_sample(bset);
  66         if (bset->n_eq == 0 && bset->n_ineq == 0)
  67                 return zero_sample(bset);
  68
  69         sample = isl_vec_alloc(bset->ctx, 2);
  70         if (!sample)
  71                 goto error;
  72         if (!bset)
  73                 return NULL;
  74         isl_int_set_si(sample->block.data[0], 1);
  75
  76         if (bset->n_eq > 0) {
  77                 isl_assert(bset->ctx, bset->n_eq == 1, goto error);
  78                 isl_assert(bset->ctx, bset->n_ineq == 0, goto error);
  79                 if (isl_int_is_one(bset->eq[0][1]))
  80                         isl_int_neg(sample->el[1], bset->eq[0][0]);
  81                 else {
  82                         isl_assert(bset->ctx, isl_int_is_negone(bset->eq[0][1]),
  83                                    goto error);
  84                         isl_int_set(sample->el[1], bset->eq[0][0]);
  85                 }
  86                 isl_basic_set_free(bset);
  87                 return sample;
  88         }
  89
  90         isl_int_init(t);
  91         if (isl_int_is_one(bset->ineq[0][1]))
  92                 isl_int_neg(sample->block.data[1], bset->ineq[0][0]);
  93         else
  94                 isl_int_set(sample->block.data[1], bset->ineq[0][0]);
  95         for (i = 1; i < bset->n_ineq; ++i) {
  96                 isl_seq_inner_product(sample->block.data,
  97                                         bset->ineq[i], 2, &t);
  98                 if (isl_int_is_neg(t))
  99                         break;
 100         }
 101         isl_int_clear(t);
 102         if (i < bset->n_ineq) {
 103                 isl_vec_free(sample);
 104                 return empty_sample(bset);
 105         }
 106
 107         isl_basic_set_free(bset);
 108         return sample;
 109 error:
 110         isl_basic_set_free(bset);
 111         isl_vec_free(sample);
 112         return NULL;
 113 }
 114
 115 /* Find a sample integer point, if any, in bset, which is known
 116  * to have equalities.  If bset contains no integer points, then
 117  * return a zero-length vector.
 118  * We simply remove the known equalities, compute a sample
 119  * in the resulting bset, using the specified recurse function,
 120  * and then transform the sample back to the original space.
 121  */
 122 static __isl_give isl_vec *sample_eq(__isl_take isl_basic_set *bset,
 123         __isl_give isl_vec *(*recurse)(__isl_take isl_basic_set *))
 124 {
 125         struct isl_mat *T;
 126         struct isl_vec *sample;
 127
 128         if (!bset)
 129                 return NULL;
 130
 131         bset = isl_basic_set_remove_equalities(bset, &T, NULL);
 132         sample = recurse(bset);
 133         if (!sample || sample->size == 0)
 134                 isl_mat_free(T);
 135         else
 136                 sample = isl_mat_vec_product(T, sample);
 137         return sample;
 138 }
 139
 140 /* Return a matrix containing the equalities of the tableau
 141  * in constraint form.  The tableau is assumed to have
 142  * an associated bset that has been kept up-to-date.
 143  */
 144 static struct isl_mat *tab_equalities(struct isl_tab *tab)
 145 {
 146         int i, j;
 147         int n_eq;
 148         struct isl_mat *eq;
 149         struct isl_basic_set *bset;
 150
 151         if (!tab)
 152                 return NULL;
 153
 154         bset = isl_tab_peek_bset(tab);
 155         isl_assert(tab->mat->ctx, bset, return NULL);
 156
 157         n_eq = tab->n_var - tab->n_col + tab->n_dead;
 158         if (tab->empty || n_eq == 0)
 159                 return isl_mat_alloc(tab->mat->ctx, 0, tab->n_var);
 160         if (n_eq == tab->n_var)
 161                 return isl_mat_identity(tab->mat->ctx, tab->n_var);
 162
 163         eq = isl_mat_alloc(tab->mat->ctx, n_eq, tab->n_var);
 164         if (!eq)
 165                 return NULL;
 166         for (i = 0, j = 0; i < tab->n_con; ++i) {
 167                 if (tab->con[i].is_row)
 168                         continue;
 169                 if (tab->con[i].index >= 0 && tab->con[i].index >= tab->n_dead)
 170                         continue;
 171                 if (i < bset->n_eq)
 172                         isl_seq_cpy(eq->row[j], bset->eq[i] + 1, tab->n_var);
 173                 else
 174                         isl_seq_cpy(eq->row[j],
 175                                     bset->ineq[i - bset->n_eq] + 1, tab->n_var);
 176                 ++j;
 177         }
 178         isl_assert(bset->ctx, j == n_eq, goto error);
 179         return eq;
 180 error:
 181         isl_mat_free(eq);
 182         return NULL;
 183 }
 184
 185 /* Compute and return an initial basis for the bounded tableau "tab".
 186  *
 187  * If the tableau is either full-dimensional or zero-dimensional,
 188  * the we simply return an identity matrix.
 189  * Otherwise, we construct a basis whose first directions correspond
 190  * to equalities.
 191  */
 192 static struct isl_mat *initial_basis(struct isl_tab *tab)
 193 {
 194         int n_eq;
 195         struct isl_mat *eq;
 196         struct isl_mat *Q;
 197
 198         tab->n_unbounded = 0;
 199         tab->n_zero = n_eq = tab->n_var - tab->n_col + tab->n_dead;
 200         if (tab->empty || n_eq == 0 || n_eq == tab->n_var)
 201                 return isl_mat_identity(tab->mat->ctx, 1 + tab->n_var);
 202
 203         eq = tab_equalities(tab);
 204         eq = isl_mat_left_hermite(eq, 0, NULL, &Q);
 205         if (!eq)
 206                 return NULL;
 207         isl_mat_free(eq);
 208
 209         Q = isl_mat_lin_to_aff(Q);
 210         return Q;
 211 }
 212
 213 /* Compute the minimum of the current ("level") basis row over "tab"
 214  * and store the result in position "level" of "min".
 215  *
 216  * This function assumes that at least one more row and at least
 217  * one more element in the constraint array are available in the tableau.
 218  */
 219 static enum isl_lp_result compute_min(isl_ctx *ctx, struct isl_tab *tab,
 220         __isl_keep isl_vec *min, int level)
 221 {
 222         return isl_tab_min(tab, tab->basis->row[1 + level],
 223                             ctx->one, &min->el[level], NULL, 0);
 224 }
 225
 226 /* Compute the maximum of the current ("level") basis row over "tab"
 227  * and store the result in position "level" of "max".
 228  *
 229  * This function assumes that at least one more row and at least
 230  * one more element in the constraint array are available in the tableau.
 231  */
 232 static enum isl_lp_result compute_max(isl_ctx *ctx, struct isl_tab *tab,
 233         __isl_keep isl_vec *max, int level)
 234 {
 235         enum isl_lp_result res;
 236         unsigned dim = tab->n_var;
 237
 238         isl_seq_neg(tab->basis->row[1 + level] + 1,
 239                     tab->basis->row[1 + level] + 1, dim);
 240         res = isl_tab_min(tab, tab->basis->row[1 + level],
 241                     ctx->one, &max->el[level], NULL, 0);
 242         isl_seq_neg(tab->basis->row[1 + level] + 1,
 243                     tab->basis->row[1 + level] + 1, dim);
 244         isl_int_neg(max->el[level], max->el[level]);
 245
 246         return res;
 247 }
 248
 249 /* Perform a greedy search for an integer point in the set represented
 250  * by "tab", given that the minimal rational value (rounded up to the
 251  * nearest integer) at "level" is smaller than the maximal rational
 252  * value (rounded down to the nearest integer).
 253  *
 254  * Return 1 if we have found an integer point (if tab->n_unbounded > 0
 255  * then we may have only found integer values for the bounded dimensions
 256  * and it is the responsibility of the caller to extend this solution
 257  * to the unbounded dimensions).
 258  * Return 0 if greedy search did not result in a solution.
 259  * Return -1 if some error occurred.
 260  *
 261  * We assign a value half-way between the minimum and the maximum
 262  * to the current dimension and check if the minimal value of the
 263  * next dimension is still smaller than (or equal) to the maximal value.
 264  * We continue this process until either
 265  * - the minimal value (rounded up) is greater than the maximal value
 266  *      (rounded down).  In this case, greedy search has failed.
 267  * - we have exhausted all bounded dimensions, meaning that we have
 268  *      found a solution.
 269  * - the sample value of the tableau is integral.
 270  * - some error has occurred.
 271  */
 272 static int greedy_search(isl_ctx *ctx, struct isl_tab *tab,
 273         __isl_keep isl_vec *min, __isl_keep isl_vec *max, int level)
 274 {
 275         struct isl_tab_undo *snap;
 276         enum isl_lp_result res;
 277
 278         snap = isl_tab_snap(tab);
 279
 280         do {
 281                 isl_int_add(tab->basis->row[1 + level][0],
 282                             min->el[level], max->el[level]);
 283                 isl_int_fdiv_q_ui(tab->basis->row[1 + level][0],
 284                             tab->basis->row[1 + level][0], 2);
 285                 isl_int_neg(tab->basis->row[1 + level][0],
 286                             tab->basis->row[1 + level][0]);
 287                 if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
 288                         return -1;
 289                 isl_int_set_si(tab->basis->row[1 + level][0], 0);
 290
 291                 if (++level >= tab->n_var - tab->n_unbounded)
 292                         return 1;
 293                 if (isl_tab_sample_is_integer(tab))
 294                         return 1;
 295
 296                 res = compute_min(ctx, tab, min, level);
 297                 if (res == isl_lp_error)
 298                         return -1;
 299                 if (res != isl_lp_ok)
 300                         isl_die(ctx, isl_error_internal,
 301                                 "expecting bounded rational solution",
 302                                 return -1);
 303                 res = compute_max(ctx, tab, max, level);
 304                 if (res == isl_lp_error)
 305                         return -1;
 306                 if (res != isl_lp_ok)
 307                         isl_die(ctx, isl_error_internal,
 308                                 "expecting bounded rational solution",
 309                                 return -1);
 310         } while (isl_int_le(min->el[level], max->el[level]));
 311
 312         if (isl_tab_rollback(tab, snap) < 0)
 313                 return -1;
 314
 315         return 0;
 316 }
 317
 318 /* Given a tableau representing a set, find and return
 319  * an integer point in the set, if there is any.
 320  *
 321  * We perform a depth first search
 322  * for an integer point, by scanning all possible values in the range
 323  * attained by a basis vector, where an initial basis may have been set
 324  * by the calling function.  Otherwise an initial basis that exploits
 325  * the equalities in the tableau is created.
 326  * tab->n_zero is currently ignored and is clobbered by this function.
 327  *
 328  * The tableau is allowed to have unbounded direction, but then
 329  * the calling function needs to set an initial basis, with the
 330  * unbounded directions last and with tab->n_unbounded set
 331  * to the number of unbounded directions.
 332  * Furthermore, the calling functions needs to add shifted copies
 333  * of all constraints involving unbounded directions to ensure
 334  * that any feasible rational value in these directions can be rounded
 335  * up to yield a feasible integer value.
 336  * In particular, let B define the given basis x' = B x
 337  * and let T be the inverse of B, i.e., X = T x'.
 338  * Let a x + c >= 0 be a constraint of the set represented by the tableau,
 339  * or a T x' + c >= 0 in terms of the given basis.  Assume that
 340  * the bounded directions have an integer value, then we can safely
 341  * round up the values for the unbounded directions if we make sure
 342  * that x' not only satisfies the original constraint, but also
 343  * the constraint "a T x' + c + s >= 0" with s the sum of all
 344  * negative values in the last n_unbounded entries of "a T".
 345  * The calling function therefore needs to add the constraint
 346  * a x + c + s >= 0.  The current function then scans the first
 347  * directions for an integer value and once those have been found,
 348  * it can compute "T ceil(B x)" to yield an integer point in the set.
 349  * Note that during the search, the first rows of B may be changed
 350  * by a basis reduction, but the last n_unbounded rows of B remain
 351  * unaltered and are also not mixed into the first rows.
 352  *
 353  * The search is implemented iteratively.  "level" identifies the current
 354  * basis vector.  "init" is true if we want the first value at the current
 355  * level and false if we want the next value.
 356  *
 357  * At the start of each level, we first check if we can find a solution
 358  * using greedy search.  If not, we continue with the exhaustive search.
 359  *
 360  * The initial basis is the identity matrix.  If the range in some direction
 361  * contains more than one integer value, we perform basis reduction based
 362  * on the value of ctx->opt->gbr
 363  *      - ISL_GBR_NEVER:        never perform basis reduction
 364  *      - ISL_GBR_ONCE:         only perform basis reduction the first
 365  *                              time such a range is encountered
 366  *      - ISL_GBR_ALWAYS:       always perform basis reduction when
 367  *                              such a range is encountered
 368  *
 369  * When ctx->opt->gbr is set to ISL_GBR_ALWAYS, then we allow the basis
 370  * reduction computation to return early.  That is, as soon as it
 371  * finds a reasonable first direction.
 372  */
 373 struct isl_vec *isl_tab_sample(struct isl_tab *tab)
 374 {
 375         unsigned dim;
 376         unsigned gbr;
 377         struct isl_ctx *ctx;
 378         struct isl_vec *sample;
 379         struct isl_vec *min;
 380         struct isl_vec *max;
 381         enum isl_lp_result res;
 382         int level;
 383         int init;
 384         int reduced;
 385         struct isl_tab_undo **snap;
 386
 387         if (!tab)
 388                 return NULL;
 389         if (tab->empty)
 390                 return isl_vec_alloc(tab->mat->ctx, 0);
 391
 392         if (!tab->basis)
 393                 tab->basis = initial_basis(tab);
 394         if (!tab->basis)
 395                 return NULL;
 396         isl_assert(tab->mat->ctx, tab->basis->n_row == tab->n_var + 1,
 397                     return NULL);
 398         isl_assert(tab->mat->ctx, tab->basis->n_col == tab->n_var + 1,
 399                     return NULL);
 400
 401         ctx = tab->mat->ctx;
 402         dim = tab->n_var;
 403         gbr = ctx->opt->gbr;
 404
 405         if (tab->n_unbounded == tab->n_var) {
 406                 sample = isl_tab_get_sample_value(tab);
 407                 sample = isl_mat_vec_product(isl_mat_copy(tab->basis), sample);
 408                 sample = isl_vec_ceil(sample);
 409                 sample = isl_mat_vec_inverse_product(isl_mat_copy(tab->basis),
 410                                                         sample);
 411                 return sample;
 412         }
 413
 414         if (isl_tab_extend_cons(tab, dim + 1) < 0)
 415                 return NULL;
 416
 417         min = isl_vec_alloc(ctx, dim);
 418         max = isl_vec_alloc(ctx, dim);
 419         snap = isl_alloc_array(ctx, struct isl_tab_undo *, dim);
 420
 421         if (!min || !max || !snap)
 422                 goto error;
 423
 424         level = 0;
 425         init = 1;
 426         reduced = 0;
 427
 428         while (level >= 0) {
 429                 if (init) {
 430                         int choice;
 431
 432                         res = compute_min(ctx, tab, min, level);
 433                         if (res == isl_lp_error)
 434                                 goto error;
 435                         if (res != isl_lp_ok)
 436                                 isl_die(ctx, isl_error_internal,
 437                                         "expecting bounded rational solution",
 438                                         goto error);
 439                         if (isl_tab_sample_is_integer(tab))
 440                                 break;
 441                         res = compute_max(ctx, tab, max, level);
 442                         if (res == isl_lp_error)
 443                                 goto error;
 444                         if (res != isl_lp_ok)
 445                                 isl_die(ctx, isl_error_internal,
 446                                         "expecting bounded rational solution",
 447                                         goto error);
 448                         if (isl_tab_sample_is_integer(tab))
 449                                 break;
 450                         choice = isl_int_lt(min->el[level], max->el[level]);
 451                         if (choice) {
 452                                 int g;
 453                                 g = greedy_search(ctx, tab, min, max, level);
 454                                 if (g < 0)
 455                                         goto error;
 456                                 if (g)
 457                                         break;
 458                         }
 459                         if (!reduced && choice &&
 460                             ctx->opt->gbr != ISL_GBR_NEVER) {
 461                                 unsigned gbr_only_first;
 462                                 if (ctx->opt->gbr == ISL_GBR_ONCE)
 463                                         ctx->opt->gbr = ISL_GBR_NEVER;
 464                                 tab->n_zero = level;
 465                                 gbr_only_first = ctx->opt->gbr_only_first;
 466                                 ctx->opt->gbr_only_first =
 467                                         ctx->opt->gbr == ISL_GBR_ALWAYS;
 468                                 tab = isl_tab_compute_reduced_basis(tab);
 469                                 ctx->opt->gbr_only_first = gbr_only_first;
 470                                 if (!tab || !tab->basis)
 471                                         goto error;
 472                                 reduced = 1;
 473                                 continue;
 474                         }
 475                         reduced = 0;
 476                         snap[level] = isl_tab_snap(tab);
 477                 } else
 478                         isl_int_add_ui(min->el[level], min->el[level], 1);
 479
 480                 if (isl_int_gt(min->el[level], max->el[level])) {
 481                         level--;
 482                         init = 0;
 483                         if (level >= 0)
 484                                 if (isl_tab_rollback(tab, snap[level]) < 0)
 485                                         goto error;
 486                         continue;
 487                 }
 488                 isl_int_neg(tab->basis->row[1 + level][0], min->el[level]);
 489                 if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
 490                         goto error;
 491                 isl_int_set_si(tab->basis->row[1 + level][0], 0);
 492                 if (level + tab->n_unbounded < dim - 1) {
 493                         ++level;
 494                         init = 1;
 495                         continue;
 496                 }
 497                 break;
 498         }
 499
 500         if (level >= 0) {
 501                 sample = isl_tab_get_sample_value(tab);
 502                 if (!sample)
 503                         goto error;
 504                 if (tab->n_unbounded && !isl_int_is_one(sample->el[0])) {
 505                         sample = isl_mat_vec_product(isl_mat_copy(tab->basis),
 506                                                      sample);
 507                         sample = isl_vec_ceil(sample);
 508                         sample = isl_mat_vec_inverse_product(
 509                                         isl_mat_copy(tab->basis), sample);
 510                 }
 511         } else
 512                 sample = isl_vec_alloc(ctx, 0);
 513
 514         ctx->opt->gbr = gbr;
 515         isl_vec_free(min);
 516         isl_vec_free(max);
 517         free(snap);
 518         return sample;
 519 error:
 520         ctx->opt->gbr = gbr;
 521         isl_vec_free(min);
 522         isl_vec_free(max);
 523         free(snap);
 524         return NULL;
 525 }
 526
 527 static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset);
 528
 529 /* Compute a sample point of the given basic set, based on the given,
 530  * non-trivial factorization.
 531  */
 532 static __isl_give isl_vec *factored_sample(__isl_take isl_basic_set *bset,
 533         __isl_take isl_factorizer *f)
 534 {
 535         int i, n;
 536         isl_vec *sample = NULL;
 537         isl_ctx *ctx;
 538         unsigned nparam;
 539         unsigned nvar;
 540
 541         ctx = isl_basic_set_get_ctx(bset);
 542         if (!ctx)
 543                 goto error;
 544
 545         nparam = isl_basic_set_dim(bset, isl_dim_param);
 546         nvar = isl_basic_set_dim(bset, isl_dim_set);
 547
 548         sample = isl_vec_alloc(ctx, 1 + isl_basic_set_total_dim(bset));
 549         if (!sample)
 550                 goto error;
 551         isl_int_set_si(sample->el[0], 1);
 552
 553         bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
 554
 555         for (i = 0, n = 0; i < f->n_group; ++i) {
 556                 isl_basic_set *bset_i;
 557                 isl_vec *sample_i;
 558
 559                 bset_i = isl_basic_set_copy(bset);
 560                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
 561                             nparam + n + f->len[i], nvar - n - f->len[i]);
 562                 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
 563                             nparam, n);
 564                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
 565                             n + f->len[i], nvar - n - f->len[i]);
 566                 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
 567
 568                 sample_i = sample_bounded(bset_i);
 569                 if (!sample_i)
 570                         goto error;
 571                 if (sample_i->size == 0) {
 572                         isl_basic_set_free(bset);
 573                         isl_factorizer_free(f);
 574                         isl_vec_free(sample);
 575                         return sample_i;
 576                 }
 577                 isl_seq_cpy(sample->el + 1 + nparam + n,
 578                             sample_i->el + 1, f->len[i]);
 579                 isl_vec_free(sample_i);
 580
 581                 n += f->len[i];
 582         }
 583
 584         f->morph = isl_morph_inverse(f->morph);
 585         sample = isl_morph_vec(isl_morph_copy(f->morph), sample);
 586
 587         isl_basic_set_free(bset);
 588         isl_factorizer_free(f);
 589         return sample;
 590 error:
 591         isl_basic_set_free(bset);
 592         isl_factorizer_free(f);
 593         isl_vec_free(sample);
 594         return NULL;
 595 }
 596
 597 /* Given a basic set that is known to be bounded, find and return
 598  * an integer point in the basic set, if there is any.
 599  *
 600  * After handling some trivial cases, we construct a tableau
 601  * and then use isl_tab_sample to find a sample, passing it
 602  * the identity matrix as initial basis.
 603  */
 604 static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset)
 605 {
 606         unsigned dim;
 607         struct isl_vec *sample;
 608         struct isl_tab *tab = NULL;
 609         isl_factorizer *f;
 610
 611         if (!bset)
 612                 return NULL;
 613
 614         if (isl_basic_set_plain_is_empty(bset))
 615                 return empty_sample(bset);
 616
 617         dim = isl_basic_set_total_dim(bset);
 618         if (dim == 0)
 619                 return zero_sample(bset);
 620         if (dim == 1)
 621                 return interval_sample(bset);
 622         if (bset->n_eq > 0)
 623                 return sample_eq(bset, sample_bounded);
 624
 625         f = isl_basic_set_factorizer(bset);
 626         if (!f)
 627                 goto error;
 628         if (f->n_group != 0)
 629                 return factored_sample(bset, f);
 630         isl_factorizer_free(f);
 631
 632         tab = isl_tab_from_basic_set(bset, 1);
 633         if (tab && tab->empty) {
 634                 isl_tab_free(tab);
 635                 ISL_F_SET(bset, ISL_BASIC_SET_EMPTY);
 636                 sample = isl_vec_alloc(isl_basic_set_get_ctx(bset), 0);
 637                 isl_basic_set_free(bset);
 638                 return sample;
 639         }
 640
 641         if (!ISL_F_ISSET(bset, ISL_BASIC_SET_NO_IMPLICIT))
 642                 if (isl_tab_detect_implicit_equalities(tab) < 0)
 643                         goto error;
 644
 645         sample = isl_tab_sample(tab);
 646         if (!sample)
 647                 goto error;
 648
 649         if (sample->size > 0) {
 650                 isl_vec_free(bset->sample);
 651                 bset->sample = isl_vec_copy(sample);
 652         }
 653
 654         isl_basic_set_free(bset);
 655         isl_tab_free(tab);
 656         return sample;
 657 error:
 658         isl_basic_set_free(bset);
 659         isl_tab_free(tab);
 660         return NULL;
 661 }
 662
 663 /* Given a basic set "bset" and a value "sample" for the first coordinates
 664  * of bset, plug in these values and drop the corresponding coordinates.
 665  *
 666  * We do this by computing the preimage of the transformation
 667  *
 668  *           [ 1 0 ]
 669  *      x =  [ s 0 ] x'
 670  *           [ 0 I ]
 671  *
 672  * where [1 s] is the sample value and I is the identity matrix of the
 673  * appropriate dimension.
 674  */
 675 static __isl_give isl_basic_set *plug_in(__isl_take isl_basic_set *bset,
 676         __isl_take isl_vec *sample)
 677 {
 678         int i;
 679         unsigned total;
 680         struct isl_mat *T;
 681
 682         if (!bset || !sample)
 683                 goto error;
 684
 685         total = isl_basic_set_total_dim(bset);
 686         T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
 687         if (!T)
 688                 goto error;
 689
 690         for (i = 0; i < sample->size; ++i) {
 691                 isl_int_set(T->row[i][0], sample->el[i]);
 692                 isl_seq_clr(T->row[i] + 1, T->n_col - 1);
 693         }
 694         for (i = 0; i < T->n_col - 1; ++i) {
 695                 isl_seq_clr(T->row[sample->size + i], T->n_col);
 696                 isl_int_set_si(T->row[sample->size + i][1 + i], 1);
 697         }
 698         isl_vec_free(sample);
 699
 700         bset = isl_basic_set_preimage(bset, T);
 701         return bset;
 702 error:
 703         isl_basic_set_free(bset);
 704         isl_vec_free(sample);
 705         return NULL;
 706 }
 707
 708 /* Given a basic set "bset", return any (possibly non-integer) point
 709  * in the basic set.
 710  */
 711 static __isl_give isl_vec *rational_sample(__isl_take isl_basic_set *bset)
 712 {
 713         struct isl_tab *tab;
 714         struct isl_vec *sample;
 715
 716         if (!bset)
 717                 return NULL;
 718
 719         tab = isl_tab_from_basic_set(bset, 0);
 720         sample = isl_tab_get_sample_value(tab);
 721         isl_tab_free(tab);
 722
 723         isl_basic_set_free(bset);
 724
 725         return sample;
 726 }
 727
 728 /* Given a linear cone "cone" and a rational point "vec",
 729  * construct a polyhedron with shifted copies of the constraints in "cone",
 730  * i.e., a polyhedron with "cone" as its recession cone, such that each
 731  * point x in this polyhedron is such that the unit box positioned at x
 732  * lies entirely inside the affine cone 'vec + cone'.
 733  * Any rational point in this polyhedron may therefore be rounded up
 734  * to yield an integer point that lies inside said affine cone.
 735  *
 736  * Denote the constraints of cone by "<a_i, x> >= 0" and the rational
 737  * point "vec" by v/d.
 738  * Let b_i = <a_i, v>.  Then the affine cone 'vec + cone' is given
 739  * by <a_i, x> - b/d >= 0.
 740  * The polyhedron <a_i, x> - ceil{b/d} >= 0 is a subset of this affine cone.
 741  * We prefer this polyhedron over the actual affine cone because it doesn't
 742  * require a scaling of the constraints.
 743  * If each of the vertices of the unit cube positioned at x lies inside
 744  * this polyhedron, then the whole unit cube at x lies inside the affine cone.
 745  * We therefore impose that x' = x + \sum e_i, for any selection of unit
 746  * vectors lies inside the polyhedron, i.e.,
 747  *
 748  *      <a_i, x'> - ceil{b/d} = <a_i, x> + sum a_i - ceil{b/d} >= 0
 749  *
 750  * The most stringent of these constraints is the one that selects
 751  * all negative a_i, so the polyhedron we are looking for has constraints
 752  *
 753  *      <a_i, x> + sum_{a_i < 0} a_i - ceil{b/d} >= 0
 754  *
 755  * Note that if cone were known to have only non-negative rays
 756  * (which can be accomplished by a unimodular transformation),
 757  * then we would only have to check the points x' = x + e_i
 758  * and we only have to add the smallest negative a_i (if any)
 759  * instead of the sum of all negative a_i.
 760  */
 761 static __isl_give isl_basic_set *shift_cone(__isl_take isl_basic_set *cone,
 762         __isl_take isl_vec *vec)
 763 {
 764         int i, j, k;
 765         unsigned total;
 766
 767         struct isl_basic_set *shift = NULL;
 768
 769         if (!cone || !vec)
 770                 goto error;
 771
 772         isl_assert(cone->ctx, cone->n_eq == 0, goto error);
 773
 774         total = isl_basic_set_total_dim(cone);
 775
 776         shift = isl_basic_set_alloc_space(isl_basic_set_get_space(cone),
 777                                         0, 0, cone->n_ineq);
 778
 779         for (i = 0; i < cone->n_ineq; ++i) {
 780                 k = isl_basic_set_alloc_inequality(shift);
 781                 if (k < 0)
 782                         goto error;
 783                 isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
 784                 isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
 785                                       &shift->ineq[k][0]);
 786                 isl_int_cdiv_q(shift->ineq[k][0],
 787                                shift->ineq[k][0], vec->el[0]);
 788                 isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
 789                 for (j = 0; j < total; ++j) {
 790                         if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
 791                                 continue;
 792                         isl_int_add(shift->ineq[k][0],
 793                                     shift->ineq[k][0], shift->ineq[k][1 + j]);
 794                 }
 795         }
 796
 797         isl_basic_set_free(cone);
 798         isl_vec_free(vec);
 799
 800         return isl_basic_set_finalize(shift);
 801 error:
 802         isl_basic_set_free(shift);
 803         isl_basic_set_free(cone);
 804         isl_vec_free(vec);
 805         return NULL;
 806 }
 807
 808 /* Given a rational point vec in a (transformed) basic set,
 809  * such that cone is the recession cone of the original basic set,
 810  * "round up" the rational point to an integer point.
 811  *
 812  * We first check if the rational point just happens to be integer.
 813  * If not, we transform the cone in the same way as the basic set,
 814  * pick a point x in this cone shifted to the rational point such that
 815  * the whole unit cube at x is also inside this affine cone.
 816  * Then we simply round up the coordinates of x and return the
 817  * resulting integer point.
 818  */
 819 static __isl_give isl_vec *round_up_in_cone(__isl_take isl_vec *vec,
 820         __isl_take isl_basic_set *cone, __isl_take isl_mat *U)
 821 {
 822         unsigned total;
 823
 824         if (!vec || !cone || !U)
 825                 goto error;
 826
 827         isl_assert(vec->ctx, vec->size != 0, goto error);
 828         if (isl_int_is_one(vec->el[0])) {
 829                 isl_mat_free(U);
 830                 isl_basic_set_free(cone);
 831                 return vec;
 832         }
 833
 834         total = isl_basic_set_total_dim(cone);
 835         cone = isl_basic_set_preimage(cone, U);
 836         cone = isl_basic_set_remove_dims(cone, isl_dim_set,
 837                                          0, total - (vec->size - 1));
 838
 839         cone = shift_cone(cone, vec);
 840
 841         vec = rational_sample(cone);
 842         vec = isl_vec_ceil(vec);
 843         return vec;
 844 error:
 845         isl_mat_free(U);
 846         isl_vec_free(vec);
 847         isl_basic_set_free(cone);
 848         return NULL;
 849 }
 850
 851 /* Concatenate two integer vectors, i.e., two vectors with denominator
 852  * (stored in element 0) equal to 1.
 853  */
 854 static __isl_give isl_vec *vec_concat(__isl_take isl_vec *vec1,
 855         __isl_take isl_vec *vec2)
 856 {
 857         struct isl_vec *vec;
 858
 859         if (!vec1 || !vec2)
 860                 goto error;
 861         isl_assert(vec1->ctx, vec1->size > 0, goto error);
 862         isl_assert(vec2->ctx, vec2->size > 0, goto error);
 863         isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);
 864         isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);
 865
 866         vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
 867         if (!vec)
 868                 goto error;
 869
 870         isl_seq_cpy(vec->el, vec1->el, vec1->size);
 871         isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
 872
 873         isl_vec_free(vec1);
 874         isl_vec_free(vec2);
 875
 876         return vec;
 877 error:
 878         isl_vec_free(vec1);
 879         isl_vec_free(vec2);
 880         return NULL;
 881 }
 882
 883 /* Give a basic set "bset" with recession cone "cone", compute and
 884  * return an integer point in bset, if any.
 885  *
 886  * If the recession cone is full-dimensional, then we know that
 887  * bset contains an infinite number of integer points and it is
 888  * fairly easy to pick one of them.
 889  * If the recession cone is not full-dimensional, then we first
 890  * transform bset such that the bounded directions appear as
 891  * the first dimensions of the transformed basic set.
 892  * We do this by using a unimodular transformation that transforms
 893  * the equalities in the recession cone to equalities on the first
 894  * dimensions.
 895  *
 896  * The transformed set is then projected onto its bounded dimensions.
 897  * Note that to compute this projection, we can simply drop all constraints
 898  * involving any of the unbounded dimensions since these constraints
 899  * cannot be combined to produce a constraint on the bounded dimensions.
 900  * To see this, assume that there is such a combination of constraints
 901  * that produces a constraint on the bounded dimensions.  This means
 902  * that some combination of the unbounded dimensions has both an upper
 903  * bound and a lower bound in terms of the bounded dimensions, but then
 904  * this combination would be a bounded direction too and would have been
 905  * transformed into a bounded dimensions.
 906  *
 907  * We then compute a sample value in the bounded dimensions.
 908  * If no such value can be found, then the original set did not contain
 909  * any integer points and we are done.
 910  * Otherwise, we plug in the value we found in the bounded dimensions,
 911  * project out these bounded dimensions and end up with a set with
 912  * a full-dimensional recession cone.
 913  * A sample point in this set is computed by "rounding up" any
 914  * rational point in the set.
 915  *
 916  * The sample points in the bounded and unbounded dimensions are
 917  * then combined into a single sample point and transformed back
 918  * to the original space.
 919  */
 920 __isl_give isl_vec *isl_basic_set_sample_with_cone(
 921         __isl_take isl_basic_set *bset, __isl_take isl_basic_set *cone)
 922 {
 923         unsigned total;
 924         unsigned cone_dim;
 925         struct isl_mat *M, *U;
 926         struct isl_vec *sample;
 927         struct isl_vec *cone_sample;
 928         struct isl_ctx *ctx;
 929         struct isl_basic_set *bounded;
 930
 931         if (!bset || !cone)
 932                 goto error;
 933
 934         ctx = isl_basic_set_get_ctx(bset);
 935         total = isl_basic_set_total_dim(cone);
 936         cone_dim = total - cone->n_eq;
 937
 938         M = isl_mat_sub_alloc6(ctx, cone->eq, 0, cone->n_eq, 1, total);
 939         M = isl_mat_left_hermite(M, 0, &U, NULL);
 940         if (!M)
 941                 goto error;
 942         isl_mat_free(M);
 943
 944         U = isl_mat_lin_to_aff(U);
 945         bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
 946
 947         bounded = isl_basic_set_copy(bset);
 948         bounded = isl_basic_set_drop_constraints_involving(bounded,
 949                                                    total - cone_dim, cone_dim);
 950         bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
 951         sample = sample_bounded(bounded);
 952         if (!sample || sample->size == 0) {
 953                 isl_basic_set_free(bset);
 954                 isl_basic_set_free(cone);
 955                 isl_mat_free(U);
 956                 return sample;
 957         }
 958         bset = plug_in(bset, isl_vec_copy(sample));
 959         cone_sample = rational_sample(bset);
 960         cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(U));
 961         sample = vec_concat(sample, cone_sample);
 962         sample = isl_mat_vec_product(U, sample);
 963         return sample;
 964 error:
 965         isl_basic_set_free(cone);
 966         isl_basic_set_free(bset);
 967         return NULL;
 968 }
 969
 970 static void vec_sum_of_neg(struct isl_vec *v, isl_int *s)
 971 {
 972         int i;
 973
 974         isl_int_set_si(*s, 0);
 975
 976         for (i = 0; i < v->size; ++i)
 977                 if (isl_int_is_neg(v->el[i]))
 978                         isl_int_add(*s, *s, v->el[i]);
 979 }
 980
 981 /* Given a tableau "tab", a tableau "tab_cone" that corresponds
 982  * to the recession cone and the inverse of a new basis U = inv(B),
 983  * with the unbounded directions in B last,
 984  * add constraints to "tab" that ensure any rational value
 985  * in the unbounded directions can be rounded up to an integer value.
 986  *
 987  * The new basis is given by x' = B x, i.e., x = U x'.
 988  * For any rational value of the last tab->n_unbounded coordinates
 989  * in the update tableau, the value that is obtained by rounding
 990  * up this value should be contained in the original tableau.
 991  * For any constraint "a x + c >= 0", we therefore need to add
 992  * a constraint "a x + c + s >= 0", with s the sum of all negative
 993  * entries in the last elements of "a U".
 994  *
 995  * Since we are not interested in the first entries of any of the "a U",
 996  * we first drop the columns of U that correpond to bounded directions.
 997  */
 998 static int tab_shift_cone(struct isl_tab *tab,
 999         struct isl_tab *tab_cone, struct isl_mat *U)
1000 {
1001         int i;
1002         isl_int v;
1003         struct isl_basic_set *bset = NULL;
1004
1005         if (tab && tab->n_unbounded == 0) {
1006                 isl_mat_free(U);
1007                 return 0;
1008         }
1009         isl_int_init(v);
1010         if (!tab || !tab_cone || !U)
1011                 goto error;
1012         bset = isl_tab_peek_bset(tab_cone);
1013         U = isl_mat_drop_cols(U, 0, tab->n_var - tab->n_unbounded);
1014         for (i = 0; i < bset->n_ineq; ++i) {
1015                 int ok;
1016                 struct isl_vec *row = NULL;
1017                 if (isl_tab_is_equality(tab_cone, tab_cone->n_eq + i))
1018                         continue;
1019                 row = isl_vec_alloc(bset->ctx, tab_cone->n_var);
1020                 if (!row)
1021                         goto error;
1022                 isl_seq_cpy(row->el, bset->ineq[i] + 1, tab_cone->n_var);
1023                 row = isl_vec_mat_product(row, isl_mat_copy(U));
1024                 if (!row)
1025                         goto error;
1026                 vec_sum_of_neg(row, &v);
1027                 isl_vec_free(row);
1028                 if (isl_int_is_zero(v))
1029                         continue;
1030                 if (isl_tab_extend_cons(tab, 1) < 0)
1031                         goto error;
1032                 isl_int_add(bset->ineq[i][0], bset->ineq[i][0], v);
1033                 ok = isl_tab_add_ineq(tab, bset->ineq[i]) >= 0;
1034                 isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], v);
1035                 if (!ok)
1036                         goto error;
1037         }
1038
1039         isl_mat_free(U);
1040         isl_int_clear(v);
1041         return 0;
1042 error:
1043         isl_mat_free(U);
1044         isl_int_clear(v);
1045         return -1;
1046 }
1047
1048 /* Compute and return an initial basis for the possibly
1049  * unbounded tableau "tab".  "tab_cone" is a tableau
1050  * for the corresponding recession cone.
1051  * Additionally, add constraints to "tab" that ensure
1052  * that any rational value for the unbounded directions
1053  * can be rounded up to an integer value.
1054  *
1055  * If the tableau is bounded, i.e., if the recession cone
1056  * is zero-dimensional, then we just use inital_basis.
1057  * Otherwise, we construct a basis whose first directions
1058  * correspond to equalities, followed by bounded directions,
1059  * i.e., equalities in the recession cone.
1060  * The remaining directions are then unbounded.
1061  */
1062 int isl_tab_set_initial_basis_with_cone(struct isl_tab *tab,
1063         struct isl_tab *tab_cone)
1064 {
1065         struct isl_mat *eq;
1066         struct isl_mat *cone_eq;
1067         struct isl_mat *U, *Q;
1068
1069         if (!tab || !tab_cone)
1070                 return -1;
1071
1072         if (tab_cone->n_col == tab_cone->n_dead) {
1073                 tab->basis = initial_basis(tab);
1074                 return tab->basis ? 0 : -1;
1075         }
1076
1077         eq = tab_equalities(tab);
1078         if (!eq)
1079                 return -1;
1080         tab->n_zero = eq->n_row;
1081         cone_eq = tab_equalities(tab_cone);
1082         eq = isl_mat_concat(eq, cone_eq);
1083         if (!eq)
1084                 return -1;
1085         tab->n_unbounded = tab->n_var - (eq->n_row - tab->n_zero);
1086         eq = isl_mat_left_hermite(eq, 0, &U, &Q);
1087         if (!eq)
1088                 return -1;
1089         isl_mat_free(eq);
1090         tab->basis = isl_mat_lin_to_aff(Q);
1091         if (tab_shift_cone(tab, tab_cone, U) < 0)
1092                 return -1;
1093         if (!tab->basis)
1094                 return -1;
1095         return 0;
1096 }
1097
1098 /* Compute and return a sample point in bset using generalized basis
1099  * reduction.  We first check if the input set has a non-trivial
1100  * recession cone.  If so, we perform some extra preprocessing in
1101  * sample_with_cone.  Otherwise, we directly perform generalized basis
1102  * reduction.
1103  */
1104 static __isl_give isl_vec *gbr_sample(__isl_take isl_basic_set *bset)
1105 {
1106         unsigned dim;
1107         struct isl_basic_set *cone;
1108
1109         dim = isl_basic_set_total_dim(bset);
1110
1111         cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
1112         if (!cone)
1113                 goto error;
1114
1115         if (cone->n_eq < dim)
1116                 return isl_basic_set_sample_with_cone(bset, cone);
1117
1118         isl_basic_set_free(cone);
1119         return sample_bounded(bset);
1120 error:
1121         isl_basic_set_free(bset);
1122         return NULL;
1123 }
1124
1125 static __isl_give isl_vec *basic_set_sample(__isl_take isl_basic_set *bset,
1126         int bounded)
1127 {
1128         struct isl_ctx *ctx;
1129         unsigned dim;
1130         if (!bset)
1131                 return NULL;
1132
1133         ctx = bset->ctx;
1134         if (isl_basic_set_plain_is_empty(bset))
1135                 return empty_sample(bset);
1136
1137         dim = isl_basic_set_n_dim(bset);
1138         if (isl_basic_set_check_no_params(bset) < 0 ||
1139             isl_basic_set_check_no_locals(bset) < 0)
1140                 goto error;
1141
1142         if (bset->sample && bset->sample->size == 1 + dim) {
1143                 int contains = isl_basic_set_contains(bset, bset->sample);
1144                 if (contains < 0)
1145                         goto error;
1146                 if (contains) {
1147                         struct isl_vec *sample = isl_vec_copy(bset->sample);
1148                         isl_basic_set_free(bset);
1149                         return sample;
1150                 }
1151         }
1152         isl_vec_free(bset->sample);
1153         bset->sample = NULL;
1154
1155         if (bset->n_eq > 0)
1156                 return sample_eq(bset, bounded ? isl_basic_set_sample_bounded
1157                                                : isl_basic_set_sample_vec);
1158         if (dim == 0)
1159                 return zero_sample(bset);
1160         if (dim == 1)
1161                 return interval_sample(bset);
1162
1163         return bounded ? sample_bounded(bset) : gbr_sample(bset);
1164 error:
1165         isl_basic_set_free(bset);
1166         return NULL;
1167 }
1168
1169 __isl_give isl_vec *isl_basic_set_sample_vec(__isl_take isl_basic_set *bset)
1170 {
1171         return basic_set_sample(bset, 0);
1172 }
1173
1174 /* Compute an integer sample in "bset", where the caller guarantees
1175  * that "bset" is bounded.
1176  */
1177 __isl_give isl_vec *isl_basic_set_sample_bounded(__isl_take isl_basic_set *bset)
1178 {
1179         return basic_set_sample(bset, 1);
1180 }
1181
1182 __isl_give isl_basic_set *isl_basic_set_from_vec(__isl_take isl_vec *vec)
1183 {
1184         int i;
1185         int k;
1186         struct isl_basic_set *bset = NULL;
1187         struct isl_ctx *ctx;
1188         unsigned dim;
1189
1190         if (!vec)
1191                 return NULL;
1192         ctx = vec->ctx;
1193         isl_assert(ctx, vec->size != 0, goto error);
1194
1195         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
1196         if (!bset)
1197                 goto error;
1198         dim = isl_basic_set_n_dim(bset);
1199         for (i = dim - 1; i >= 0; --i) {
1200                 k = isl_basic_set_alloc_equality(bset);
1201                 if (k < 0)
1202                         goto error;
1203                 isl_seq_clr(bset->eq[k], 1 + dim);
1204                 isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
1205                 isl_int_set(bset->eq[k][1 + i], vec->el[0]);
1206         }
1207         bset->sample = vec;
1208
1209         return bset;
1210 error:
1211         isl_basic_set_free(bset);
1212         isl_vec_free(vec);
1213         return NULL;
1214 }
1215
1216 __isl_give isl_basic_map *isl_basic_map_sample(__isl_take isl_basic_map *bmap)
1217 {
1218         struct isl_basic_set *bset;
1219         struct isl_vec *sample_vec;
1220
1221         bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1222         sample_vec = isl_basic_set_sample_vec(bset);
1223         if (!sample_vec)
1224                 goto error;
1225         if (sample_vec->size == 0) {
1226                 isl_vec_free(sample_vec);
1227                 return isl_basic_map_set_to_empty(bmap);
1228         }
1229         isl_vec_free(bmap->sample);
1230         bmap->sample = isl_vec_copy(sample_vec);
1231         bset = isl_basic_set_from_vec(sample_vec);
1232         return isl_basic_map_overlying_set(bset, bmap);
1233 error:
1234         isl_basic_map_free(bmap);
1235         return NULL;
1236 }
1237
1238 __isl_give isl_basic_set *isl_basic_set_sample(__isl_take isl_basic_set *bset)
1239 {
1240         return isl_basic_map_sample(bset);
1241 }
1242
1243 __isl_give isl_basic_map *isl_map_sample(__isl_take isl_map *map)
1244 {
1245         int i;
1246         isl_basic_map *sample = NULL;
1247
1248         if (!map)
1249                 goto error;
1250
1251         for (i = 0; i < map->n; ++i) {
1252                 sample = isl_basic_map_sample(isl_basic_map_copy(map->p[i]));
1253                 if (!sample)
1254                         goto error;
1255                 if (!ISL_F_ISSET(sample, ISL_BASIC_MAP_EMPTY))
1256                         break;
1257                 isl_basic_map_free(sample);
1258         }
1259         if (i == map->n)
1260                 sample = isl_basic_map_empty(isl_map_get_space(map));
1261         isl_map_free(map);
1262         return sample;
1263 error:
1264         isl_map_free(map);
1265         return NULL;
1266 }
1267
1268 __isl_give isl_basic_set *isl_set_sample(__isl_take isl_set *set)
1269 {
1270         return bset_from_bmap(isl_map_sample(set_to_map(set)));
1271 }
1272
1273 __isl_give isl_point *isl_basic_set_sample_point(__isl_take isl_basic_set *bset)
1274 {
1275         isl_vec *vec;
1276         isl_space *space;
1277
1278         space = isl_basic_set_get_space(bset);
1279         bset = isl_basic_set_underlying_set(bset);
1280         vec = isl_basic_set_sample_vec(bset);
1281
1282         return isl_point_alloc(space, vec);
1283 }
1284
1285 __isl_give isl_point *isl_set_sample_point(__isl_take isl_set *set)
1286 {
1287         int i;
1288         isl_point *pnt;
1289
1290         if (!set)
1291                 return NULL;
1292
1293         for (i = 0; i < set->n; ++i) {
1294                 pnt = isl_basic_set_sample_point(isl_basic_set_copy(set->p[i]));
1295                 if (!pnt)
1296                         goto error;
1297                 if (!isl_point_is_void(pnt))
1298                         break;
1299                 isl_point_free(pnt);
1300         }
1301         if (i == set->n)
1302                 pnt = isl_point_void(isl_set_get_space(set));
1303
1304         isl_set_free(set);
1305         return pnt;
1306 error:
1307         isl_set_free(set);
1308         return NULL;
1309 }