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