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