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 __isl_give 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 /* Internal data for factored_sample.
 535  * "sample" collects the sample and may get reset to a zero-length vector
 536  * signaling the absence of a sample vector.
 537  * "pos" is the position of the contribution of the next factor.
 538  */
 539 struct isl_factored_sample_data {
 540         isl_vec *sample;
 541         int pos;
 542 };
 543
 544 /* isl_factorizer_every_factor_basic_set callback that extends
 545  * the sample in data->sample with the contribution
 546  * of the factor "bset".
 547  * If "bset" turns out to be empty, then the product is empty too and
 548  * no further factors need to be considered.
 549  */
 550 static isl_bool factor_sample(__isl_keep isl_basic_set *bset, void *user)
 551 {
 552         struct isl_factored_sample_data *data = user;
 553         isl_vec *sample;
 554         isl_size n;
 555
 556         n = isl_basic_set_dim(bset, isl_dim_set);
 557         if (n < 0)
 558                 return isl_bool_error;
 559
 560         sample = sample_bounded(isl_basic_set_copy(bset));
 561         if (!sample)
 562                 return isl_bool_error;
 563         if (sample->size == 0) {
 564                 isl_vec_free(data->sample);
 565                 data->sample = sample;
 566                 return isl_bool_false;
 567         }
 568         isl_seq_cpy(data->sample->el + data->pos, sample->el + 1, n);
 569         isl_vec_free(sample);
 570         data->pos += n;
 571
 572         return isl_bool_true;
 573 }
 574
 575 /* Compute a sample point of the given basic set, based on the given,
 576  * non-trivial factorization.
 577  */
 578 static __isl_give isl_vec *factored_sample(__isl_take isl_basic_set *bset,
 579         __isl_take isl_factorizer *f)
 580 {
 581         struct isl_factored_sample_data data = { NULL };
 582         isl_ctx *ctx;
 583         isl_size total;
 584         isl_bool every;
 585
 586         ctx = isl_basic_set_get_ctx(bset);
 587         total = isl_basic_set_dim(bset, isl_dim_all);
 588         if (!ctx || total < 0)
 589                 goto error;
 590
 591         data.sample = isl_vec_alloc(ctx, 1 + total);
 592         if (!data.sample)
 593                 goto error;
 594         isl_int_set_si(data.sample->el[0], 1);
 595         data.pos = 1;
 596
 597         every = isl_factorizer_every_factor_basic_set(f, &factor_sample, &data);
 598         if (every < 0) {
 599                 data.sample = isl_vec_free(data.sample);
 600         } else if (every) {
 601                 isl_morph *morph;
 602
 603                 morph = isl_morph_inverse(isl_morph_copy(f->morph));
 604                 data.sample = isl_morph_vec(morph, data.sample);
 605         }
 606
 607         isl_basic_set_free(bset);
 608         isl_factorizer_free(f);
 609         return data.sample;
 610 error:
 611         isl_basic_set_free(bset);
 612         isl_factorizer_free(f);
 613         isl_vec_free(data.sample);
 614         return NULL;
 615 }
 616
 617 /* Given a basic set that is known to be bounded, find and return
 618  * an integer point in the basic set, if there is any.
 619  *
 620  * After handling some trivial cases, we construct a tableau
 621  * and then use isl_tab_sample to find a sample, passing it
 622  * the identity matrix as initial basis.
 623  */
 624 static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset)
 625 {
 626         isl_size dim;
 627         struct isl_vec *sample;
 628         struct isl_tab *tab = NULL;
 629         isl_factorizer *f;
 630
 631         if (!bset)
 632                 return NULL;
 633
 634         if (isl_basic_set_plain_is_empty(bset))
 635                 return empty_sample(bset);
 636
 637         dim = isl_basic_set_dim(bset, isl_dim_all);
 638         if (dim < 0)
 639                 bset = isl_basic_set_free(bset);
 640         if (dim == 0)
 641                 return zero_sample(bset);
 642         if (dim == 1)
 643                 return interval_sample(bset);
 644         if (bset->n_eq > 0)
 645                 return sample_eq(bset, sample_bounded);
 646
 647         f = isl_basic_set_factorizer(bset);
 648         if (!f)
 649                 goto error;
 650         if (f->n_group != 0)
 651                 return factored_sample(bset, f);
 652         isl_factorizer_free(f);
 653
 654         tab = isl_tab_from_basic_set(bset, 1);
 655         if (tab && tab->empty) {
 656                 isl_tab_free(tab);
 657                 ISL_F_SET(bset, ISL_BASIC_SET_EMPTY);
 658                 sample = isl_vec_alloc(isl_basic_set_get_ctx(bset), 0);
 659                 isl_basic_set_free(bset);
 660                 return sample;
 661         }
 662
 663         if (!ISL_F_ISSET(bset, ISL_BASIC_SET_NO_IMPLICIT))
 664                 if (isl_tab_detect_implicit_equalities(tab) < 0)
 665                         goto error;
 666
 667         sample = isl_tab_sample(tab);
 668         if (!sample)
 669                 goto error;
 670
 671         if (sample->size > 0) {
 672                 isl_vec_free(bset->sample);
 673                 bset->sample = isl_vec_copy(sample);
 674         }
 675
 676         isl_basic_set_free(bset);
 677         isl_tab_free(tab);
 678         return sample;
 679 error:
 680         isl_basic_set_free(bset);
 681         isl_tab_free(tab);
 682         return NULL;
 683 }
 684
 685 /* Given a basic set "bset" and a value "sample" for the first coordinates
 686  * of bset, plug in these values and drop the corresponding coordinates.
 687  *
 688  * We do this by computing the preimage of the transformation
 689  *
 690  *           [ 1 0 ]
 691  *      x =  [ s 0 ] x'
 692  *           [ 0 I ]
 693  *
 694  * where [1 s] is the sample value and I is the identity matrix of the
 695  * appropriate dimension.
 696  */
 697 static __isl_give isl_basic_set *plug_in(__isl_take isl_basic_set *bset,
 698         __isl_take isl_vec *sample)
 699 {
 700         int i;
 701         isl_size total;
 702         struct isl_mat *T;
 703
 704         total = isl_basic_set_dim(bset, isl_dim_all);
 705         if (total < 0 || !sample)
 706                 goto error;
 707
 708         T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
 709         if (!T)
 710                 goto error;
 711
 712         for (i = 0; i < sample->size; ++i) {
 713                 isl_int_set(T->row[i][0], sample->el[i]);
 714                 isl_seq_clr(T->row[i] + 1, T->n_col - 1);
 715         }
 716         for (i = 0; i < T->n_col - 1; ++i) {
 717                 isl_seq_clr(T->row[sample->size + i], T->n_col);
 718                 isl_int_set_si(T->row[sample->size + i][1 + i], 1);
 719         }
 720         isl_vec_free(sample);
 721
 722         bset = isl_basic_set_preimage(bset, T);
 723         return bset;
 724 error:
 725         isl_basic_set_free(bset);
 726         isl_vec_free(sample);
 727         return NULL;
 728 }
 729
 730 /* Given a basic set "bset", return any (possibly non-integer) point
 731  * in the basic set.
 732  */
 733 static __isl_give isl_vec *rational_sample(__isl_take isl_basic_set *bset)
 734 {
 735         struct isl_tab *tab;
 736         struct isl_vec *sample;
 737
 738         if (!bset)
 739                 return NULL;
 740
 741         tab = isl_tab_from_basic_set(bset, 0);
 742         sample = isl_tab_get_sample_value(tab);
 743         isl_tab_free(tab);
 744
 745         isl_basic_set_free(bset);
 746
 747         return sample;
 748 }
 749
 750 /* Given a linear cone "cone" and a rational point "vec",
 751  * construct a polyhedron with shifted copies of the constraints in "cone",
 752  * i.e., a polyhedron with "cone" as its recession cone, such that each
 753  * point x in this polyhedron is such that the unit box positioned at x
 754  * lies entirely inside the affine cone 'vec + cone'.
 755  * Any rational point in this polyhedron may therefore be rounded up
 756  * to yield an integer point that lies inside said affine cone.
 757  *
 758  * Denote the constraints of cone by "<a_i, x> >= 0" and the rational
 759  * point "vec" by v/d.
 760  * Let b_i = <a_i, v>.  Then the affine cone 'vec + cone' is given
 761  * by <a_i, x> - b/d >= 0.
 762  * The polyhedron <a_i, x> - ceil{b/d} >= 0 is a subset of this affine cone.
 763  * We prefer this polyhedron over the actual affine cone because it doesn't
 764  * require a scaling of the constraints.
 765  * If each of the vertices of the unit cube positioned at x lies inside
 766  * this polyhedron, then the whole unit cube at x lies inside the affine cone.
 767  * We therefore impose that x' = x + \sum e_i, for any selection of unit
 768  * vectors lies inside the polyhedron, i.e.,
 769  *
 770  *      <a_i, x'> - ceil{b/d} = <a_i, x> + sum a_i - ceil{b/d} >= 0
 771  *
 772  * The most stringent of these constraints is the one that selects
 773  * all negative a_i, so the polyhedron we are looking for has constraints
 774  *
 775  *      <a_i, x> + sum_{a_i < 0} a_i - ceil{b/d} >= 0
 776  *
 777  * Note that if cone were known to have only non-negative rays
 778  * (which can be accomplished by a unimodular transformation),
 779  * then we would only have to check the points x' = x + e_i
 780  * and we only have to add the smallest negative a_i (if any)
 781  * instead of the sum of all negative a_i.
 782  */
 783 static __isl_give isl_basic_set *shift_cone(__isl_take isl_basic_set *cone,
 784         __isl_take isl_vec *vec)
 785 {
 786         int i, j, k;
 787         isl_size total;
 788
 789         struct isl_basic_set *shift = NULL;
 790
 791         total = isl_basic_set_dim(cone, isl_dim_all);
 792         if (total < 0 || !vec)
 793                 goto error;
 794
 795         isl_assert(cone->ctx, cone->n_eq == 0, goto error);
 796
 797         shift = isl_basic_set_alloc_space(isl_basic_set_get_space(cone),
 798                                         0, 0, cone->n_ineq);
 799
 800         for (i = 0; i < cone->n_ineq; ++i) {
 801                 k = isl_basic_set_alloc_inequality(shift);
 802                 if (k < 0)
 803                         goto error;
 804                 isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
 805                 isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
 806                                       &shift->ineq[k][0]);
 807                 isl_int_cdiv_q(shift->ineq[k][0],
 808                                shift->ineq[k][0], vec->el[0]);
 809                 isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
 810                 for (j = 0; j < total; ++j) {
 811                         if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
 812                                 continue;
 813                         isl_int_add(shift->ineq[k][0],
 814                                     shift->ineq[k][0], shift->ineq[k][1 + j]);
 815                 }
 816         }
 817
 818         isl_basic_set_free(cone);
 819         isl_vec_free(vec);
 820
 821         return isl_basic_set_finalize(shift);
 822 error:
 823         isl_basic_set_free(shift);
 824         isl_basic_set_free(cone);
 825         isl_vec_free(vec);
 826         return NULL;
 827 }
 828
 829 /* Given a rational point vec in a (transformed) basic set,
 830  * such that cone is the recession cone of the original basic set,
 831  * "round up" the rational point to an integer point.
 832  *
 833  * We first check if the rational point just happens to be integer.
 834  * If not, we transform the cone in the same way as the basic set,
 835  * pick a point x in this cone shifted to the rational point such that
 836  * the whole unit cube at x is also inside this affine cone.
 837  * Then we simply round up the coordinates of x and return the
 838  * resulting integer point.
 839  */
 840 static __isl_give isl_vec *round_up_in_cone(__isl_take isl_vec *vec,
 841         __isl_take isl_basic_set *cone, __isl_take isl_mat *U)
 842 {
 843         isl_size total;
 844
 845         if (!vec || !cone || !U)
 846                 goto error;
 847
 848         isl_assert(vec->ctx, vec->size != 0, goto error);
 849         if (isl_int_is_one(vec->el[0])) {
 850                 isl_mat_free(U);
 851                 isl_basic_set_free(cone);
 852                 return vec;
 853         }
 854
 855         total = isl_basic_set_dim(cone, isl_dim_all);
 856         if (total < 0)
 857                 goto error;
 858         cone = isl_basic_set_preimage(cone, U);
 859         cone = isl_basic_set_remove_dims(cone, isl_dim_set,
 860                                          0, total - (vec->size - 1));
 861
 862         cone = shift_cone(cone, vec);
 863
 864         vec = rational_sample(cone);
 865         vec = isl_vec_ceil(vec);
 866         return vec;
 867 error:
 868         isl_mat_free(U);
 869         isl_vec_free(vec);
 870         isl_basic_set_free(cone);
 871         return NULL;
 872 }
 873
 874 /* Concatenate two integer vectors, i.e., two vectors with denominator
 875  * (stored in element 0) equal to 1.
 876  */
 877 static __isl_give isl_vec *vec_concat(__isl_take isl_vec *vec1,
 878         __isl_take isl_vec *vec2)
 879 {
 880         struct isl_vec *vec;
 881
 882         if (!vec1 || !vec2)
 883                 goto error;
 884         isl_assert(vec1->ctx, vec1->size > 0, goto error);
 885         isl_assert(vec2->ctx, vec2->size > 0, goto error);
 886         isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);
 887         isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);
 888
 889         vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
 890         if (!vec)
 891                 goto error;
 892
 893         isl_seq_cpy(vec->el, vec1->el, vec1->size);
 894         isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
 895
 896         isl_vec_free(vec1);
 897         isl_vec_free(vec2);
 898
 899         return vec;
 900 error:
 901         isl_vec_free(vec1);
 902         isl_vec_free(vec2);
 903         return NULL;
 904 }
 905
 906 /* Give a basic set "bset" with recession cone "cone", compute and
 907  * return an integer point in bset, if any.
 908  *
 909  * If the recession cone is full-dimensional, then we know that
 910  * bset contains an infinite number of integer points and it is
 911  * fairly easy to pick one of them.
 912  * If the recession cone is not full-dimensional, then we first
 913  * transform bset such that the bounded directions appear as
 914  * the first dimensions of the transformed basic set.
 915  * We do this by using a unimodular transformation that transforms
 916  * the equalities in the recession cone to equalities on the first
 917  * dimensions.
 918  *
 919  * The transformed set is then projected onto its bounded dimensions.
 920  * Note that to compute this projection, we can simply drop all constraints
 921  * involving any of the unbounded dimensions since these constraints
 922  * cannot be combined to produce a constraint on the bounded dimensions.
 923  * To see this, assume that there is such a combination of constraints
 924  * that produces a constraint on the bounded dimensions.  This means
 925  * that some combination of the unbounded dimensions has both an upper
 926  * bound and a lower bound in terms of the bounded dimensions, but then
 927  * this combination would be a bounded direction too and would have been
 928  * transformed into a bounded dimensions.
 929  *
 930  * We then compute a sample value in the bounded dimensions.
 931  * If no such value can be found, then the original set did not contain
 932  * any integer points and we are done.
 933  * Otherwise, we plug in the value we found in the bounded dimensions,
 934  * project out these bounded dimensions and end up with a set with
 935  * a full-dimensional recession cone.
 936  * A sample point in this set is computed by "rounding up" any
 937  * rational point in the set.
 938  *
 939  * The sample points in the bounded and unbounded dimensions are
 940  * then combined into a single sample point and transformed back
 941  * to the original space.
 942  */
 943 __isl_give isl_vec *isl_basic_set_sample_with_cone(
 944         __isl_take isl_basic_set *bset, __isl_take isl_basic_set *cone)
 945 {
 946         isl_size total;
 947         unsigned cone_dim;
 948         struct isl_mat *M, *U;
 949         struct isl_vec *sample;
 950         struct isl_vec *cone_sample;
 951         struct isl_ctx *ctx;
 952         struct isl_basic_set *bounded;
 953
 954         total = isl_basic_set_dim(cone, isl_dim_all);
 955         if (!bset || total < 0)
 956                 goto error;
 957
 958         ctx = isl_basic_set_get_ctx(bset);
 959         cone_dim = total - cone->n_eq;
 960
 961         M = isl_mat_sub_alloc6(ctx, cone->eq, 0, cone->n_eq, 1, total);
 962         M = isl_mat_left_hermite(M, 0, &U, NULL);
 963         if (!M)
 964                 goto error;
 965         isl_mat_free(M);
 966
 967         U = isl_mat_lin_to_aff(U);
 968         bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
 969
 970         bounded = isl_basic_set_copy(bset);
 971         bounded = isl_basic_set_drop_constraints_involving(bounded,
 972                                                    total - cone_dim, cone_dim);
 973         bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
 974         sample = sample_bounded(bounded);
 975         if (!sample || sample->size == 0) {
 976                 isl_basic_set_free(bset);
 977                 isl_basic_set_free(cone);
 978                 isl_mat_free(U);
 979                 return sample;
 980         }
 981         bset = plug_in(bset, isl_vec_copy(sample));
 982         cone_sample = rational_sample(bset);
 983         cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(U));
 984         sample = vec_concat(sample, cone_sample);
 985         sample = isl_mat_vec_product(U, sample);
 986         return sample;
 987 error:
 988         isl_basic_set_free(cone);
 989         isl_basic_set_free(bset);
 990         return NULL;
 991 }
 992
 993 static void vec_sum_of_neg(__isl_keep isl_vec *v, isl_int *s)
 994 {
 995         int i;
 996
 997         isl_int_set_si(*s, 0);
 998
 999         for (i = 0; i < v->size; ++i)
1000                 if (isl_int_is_neg(v->el[i]))
1001                         isl_int_add(*s, *s, v->el[i]);
1002 }
1003
1004 /* Given a tableau "tab", a tableau "tab_cone" that corresponds
1005  * to the recession cone and the inverse of a new basis U = inv(B),
1006  * with the unbounded directions in B last,
1007  * add constraints to "tab" that ensure any rational value
1008  * in the unbounded directions can be rounded up to an integer value.
1009  *
1010  * The new basis is given by x' = B x, i.e., x = U x'.
1011  * For any rational value of the last tab->n_unbounded coordinates
1012  * in the update tableau, the value that is obtained by rounding
1013  * up this value should be contained in the original tableau.
1014  * For any constraint "a x + c >= 0", we therefore need to add
1015  * a constraint "a x + c + s >= 0", with s the sum of all negative
1016  * entries in the last elements of "a U".
1017  *
1018  * Since we are not interested in the first entries of any of the "a U",
1019  * we first drop the columns of U that correpond to bounded directions.
1020  */
1021 static int tab_shift_cone(struct isl_tab *tab,
1022         struct isl_tab *tab_cone, struct isl_mat *U)
1023 {
1024         int i;
1025         isl_int v;
1026         struct isl_basic_set *bset = NULL;
1027
1028         if (tab && tab->n_unbounded == 0) {
1029                 isl_mat_free(U);
1030                 return 0;
1031         }
1032         isl_int_init(v);
1033         if (!tab || !tab_cone || !U)
1034                 goto error;
1035         bset = isl_tab_peek_bset(tab_cone);
1036         U = isl_mat_drop_cols(U, 0, tab->n_var - tab->n_unbounded);
1037         for (i = 0; i < bset->n_ineq; ++i) {
1038                 int ok;
1039                 struct isl_vec *row = NULL;
1040                 if (isl_tab_is_equality(tab_cone, tab_cone->n_eq + i))
1041                         continue;
1042                 row = isl_vec_alloc(bset->ctx, tab_cone->n_var);
1043                 if (!row)
1044                         goto error;
1045                 isl_seq_cpy(row->el, bset->ineq[i] + 1, tab_cone->n_var);
1046                 row = isl_vec_mat_product(row, isl_mat_copy(U));
1047                 if (!row)
1048                         goto error;
1049                 vec_sum_of_neg(row, &v);
1050                 isl_vec_free(row);
1051                 if (isl_int_is_zero(v))
1052                         continue;
1053                 if (isl_tab_extend_cons(tab, 1) < 0)
1054                         goto error;
1055                 isl_int_add(bset->ineq[i][0], bset->ineq[i][0], v);
1056                 ok = isl_tab_add_ineq(tab, bset->ineq[i]) >= 0;
1057                 isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], v);
1058                 if (!ok)
1059                         goto error;
1060         }
1061
1062         isl_mat_free(U);
1063         isl_int_clear(v);
1064         return 0;
1065 error:
1066         isl_mat_free(U);
1067         isl_int_clear(v);
1068         return -1;
1069 }
1070
1071 /* Compute and return an initial basis for the possibly
1072  * unbounded tableau "tab".  "tab_cone" is a tableau
1073  * for the corresponding recession cone.
1074  * Additionally, add constraints to "tab" that ensure
1075  * that any rational value for the unbounded directions
1076  * can be rounded up to an integer value.
1077  *
1078  * If the tableau is bounded, i.e., if the recession cone
1079  * is zero-dimensional, then we just use inital_basis.
1080  * Otherwise, we construct a basis whose first directions
1081  * correspond to equalities, followed by bounded directions,
1082  * i.e., equalities in the recession cone.
1083  * The remaining directions are then unbounded.
1084  */
1085 int isl_tab_set_initial_basis_with_cone(struct isl_tab *tab,
1086         struct isl_tab *tab_cone)
1087 {
1088         struct isl_mat *eq;
1089         struct isl_mat *cone_eq;
1090         struct isl_mat *U, *Q;
1091
1092         if (!tab || !tab_cone)
1093                 return -1;
1094
1095         if (tab_cone->n_col == tab_cone->n_dead) {
1096                 tab->basis = initial_basis(tab);
1097                 return tab->basis ? 0 : -1;
1098         }
1099
1100         eq = tab_equalities(tab);
1101         if (!eq)
1102                 return -1;
1103         tab->n_zero = eq->n_row;
1104         cone_eq = tab_equalities(tab_cone);
1105         eq = isl_mat_concat(eq, cone_eq);
1106         if (!eq)
1107                 return -1;
1108         tab->n_unbounded = tab->n_var - (eq->n_row - tab->n_zero);
1109         eq = isl_mat_left_hermite(eq, 0, &U, &Q);
1110         if (!eq)
1111                 return -1;
1112         isl_mat_free(eq);
1113         tab->basis = isl_mat_lin_to_aff(Q);
1114         if (tab_shift_cone(tab, tab_cone, U) < 0)
1115                 return -1;
1116         if (!tab->basis)
1117                 return -1;
1118         return 0;
1119 }
1120
1121 /* Compute and return a sample point in bset using generalized basis
1122  * reduction.  We first check if the input set has a non-trivial
1123  * recession cone.  If so, we perform some extra preprocessing in
1124  * sample_with_cone.  Otherwise, we directly perform generalized basis
1125  * reduction.
1126  */
1127 static __isl_give isl_vec *gbr_sample(__isl_take isl_basic_set *bset)
1128 {
1129         isl_size dim;
1130         struct isl_basic_set *cone;
1131
1132         dim = isl_basic_set_dim(bset, isl_dim_all);
1133         if (dim < 0)
1134                 goto error;
1135
1136         cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
1137         if (!cone)
1138                 goto error;
1139
1140         if (cone->n_eq < dim)
1141                 return isl_basic_set_sample_with_cone(bset, cone);
1142
1143         isl_basic_set_free(cone);
1144         return sample_bounded(bset);
1145 error:
1146         isl_basic_set_free(bset);
1147         return NULL;
1148 }
1149
1150 static __isl_give isl_vec *basic_set_sample(__isl_take isl_basic_set *bset,
1151         int bounded)
1152 {
1153         struct isl_ctx *ctx;
1154         isl_size dim;
1155         if (!bset)
1156                 return NULL;
1157
1158         ctx = bset->ctx;
1159         if (isl_basic_set_plain_is_empty(bset))
1160                 return empty_sample(bset);
1161
1162         dim = isl_basic_set_dim(bset, isl_dim_set);
1163         if (dim < 0 ||
1164             isl_basic_set_check_no_params(bset) < 0 ||
1165             isl_basic_set_check_no_locals(bset) < 0)
1166                 goto error;
1167
1168         if (bset->sample && bset->sample->size == 1 + dim) {
1169                 int contains = isl_basic_set_contains(bset, bset->sample);
1170                 if (contains < 0)
1171                         goto error;
1172                 if (contains) {
1173                         struct isl_vec *sample = isl_vec_copy(bset->sample);
1174                         isl_basic_set_free(bset);
1175                         return sample;
1176                 }
1177         }
1178         isl_vec_free(bset->sample);
1179         bset->sample = NULL;
1180
1181         if (bset->n_eq > 0)
1182                 return sample_eq(bset, bounded ? isl_basic_set_sample_bounded
1183                                                : isl_basic_set_sample_vec);
1184         if (dim == 0)
1185                 return zero_sample(bset);
1186         if (dim == 1)
1187                 return interval_sample(bset);
1188
1189         return bounded ? sample_bounded(bset) : gbr_sample(bset);
1190 error:
1191         isl_basic_set_free(bset);
1192         return NULL;
1193 }
1194
1195 __isl_give isl_vec *isl_basic_set_sample_vec(__isl_take isl_basic_set *bset)
1196 {
1197         return basic_set_sample(bset, 0);
1198 }
1199
1200 /* Compute an integer sample in "bset", where the caller guarantees
1201  * that "bset" is bounded.
1202  */
1203 __isl_give isl_vec *isl_basic_set_sample_bounded(__isl_take isl_basic_set *bset)
1204 {
1205         return basic_set_sample(bset, 1);
1206 }
1207
1208 __isl_give isl_basic_set *isl_basic_set_from_vec(__isl_take isl_vec *vec)
1209 {
1210         int i;
1211         int k;
1212         struct isl_basic_set *bset = NULL;
1213         struct isl_ctx *ctx;
1214         isl_size dim;
1215
1216         if (!vec)
1217                 return NULL;
1218         ctx = vec->ctx;
1219         isl_assert(ctx, vec->size != 0, goto error);
1220
1221         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
1222         dim = isl_basic_set_dim(bset, isl_dim_set);
1223         if (dim < 0)
1224                 goto error;
1225         for (i = dim - 1; i >= 0; --i) {
1226                 k = isl_basic_set_alloc_equality(bset);
1227                 if (k < 0)
1228                         goto error;
1229                 isl_seq_clr(bset->eq[k], 1 + dim);
1230                 isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
1231                 isl_int_set(bset->eq[k][1 + i], vec->el[0]);
1232         }
1233         bset->sample = vec;
1234
1235         return bset;
1236 error:
1237         isl_basic_set_free(bset);
1238         isl_vec_free(vec);
1239         return NULL;
1240 }
1241
1242 __isl_give isl_basic_map *isl_basic_map_sample(__isl_take isl_basic_map *bmap)
1243 {
1244         struct isl_basic_set *bset;
1245         struct isl_vec *sample_vec;
1246
1247         bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1248         sample_vec = isl_basic_set_sample_vec(bset);
1249         if (!sample_vec)
1250                 goto error;
1251         if (sample_vec->size == 0) {
1252                 isl_vec_free(sample_vec);
1253                 return isl_basic_map_set_to_empty(bmap);
1254         }
1255         isl_vec_free(bmap->sample);
1256         bmap->sample = isl_vec_copy(sample_vec);
1257         bset = isl_basic_set_from_vec(sample_vec);
1258         return isl_basic_map_overlying_set(bset, bmap);
1259 error:
1260         isl_basic_map_free(bmap);
1261         return NULL;
1262 }
1263
1264 __isl_give isl_basic_set *isl_basic_set_sample(__isl_take isl_basic_set *bset)
1265 {
1266         return isl_basic_map_sample(bset);
1267 }
1268
1269 __isl_give isl_basic_map *isl_map_sample(__isl_take isl_map *map)
1270 {
1271         int i;
1272         isl_basic_map *sample = NULL;
1273
1274         if (!map)
1275                 goto error;
1276
1277         for (i = 0; i < map->n; ++i) {
1278                 sample = isl_basic_map_sample(isl_basic_map_copy(map->p[i]));
1279                 if (!sample)
1280                         goto error;
1281                 if (!ISL_F_ISSET(sample, ISL_BASIC_MAP_EMPTY))
1282                         break;
1283                 isl_basic_map_free(sample);
1284         }
1285         if (i == map->n)
1286                 sample = isl_basic_map_empty(isl_map_get_space(map));
1287         isl_map_free(map);
1288         return sample;
1289 error:
1290         isl_map_free(map);
1291         return NULL;
1292 }
1293
1294 __isl_give isl_basic_set *isl_set_sample(__isl_take isl_set *set)
1295 {
1296         return bset_from_bmap(isl_map_sample(set_to_map(set)));
1297 }
1298
1299 __isl_give isl_point *isl_basic_set_sample_point(__isl_take isl_basic_set *bset)
1300 {
1301         isl_vec *vec;
1302         isl_space *space;
1303
1304         space = isl_basic_set_get_space(bset);
1305         bset = isl_basic_set_underlying_set(bset);
1306         vec = isl_basic_set_sample_vec(bset);
1307
1308         return isl_point_alloc(space, vec);
1309 }
1310
1311 __isl_give isl_point *isl_set_sample_point(__isl_take isl_set *set)
1312 {
1313         int i;
1314         isl_point *pnt;
1315
1316         if (!set)
1317                 return NULL;
1318
1319         for (i = 0; i < set->n; ++i) {
1320                 pnt = isl_basic_set_sample_point(isl_basic_set_copy(set->p[i]));
1321                 if (!pnt)
1322                         goto error;
1323                 if (!isl_point_is_void(pnt))
1324                         break;
1325                 isl_point_free(pnt);
1326         }
1327         if (i == set->n)
1328                 pnt = isl_point_void(isl_set_get_space(set));
1329
1330         isl_set_free(set);
1331         return pnt;
1332 error:
1333         isl_set_free(set);
1334         return NULL;
1335 }