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