optionally use imath for arbitrary precision integers
[isl.git] / basis_reduction_tab.c
blob112d423dbe0987f3a20d0b72a4d2ca22fe774671
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  */
10 #include <assert.h>
11 #include <isl_map_private.h>
12 #include <isl_seq.h>
13 #include "isl_tab.h"
14 #include <isl_int.h>
15 #include <isl_config.h>
17 struct tab_lp {
18         struct isl_ctx  *ctx;
19         struct isl_vec  *row;
20         struct isl_tab  *tab;
21         struct isl_tab_undo     **stack;
22         isl_int         *obj;
23         isl_int          opt;
24         isl_int          opt_denom;
25         isl_int          tmp;
26         isl_int          tmp2;
27         int              neq;
28         unsigned         dim;
29         /* number of constraints in initial product tableau */
30         int              con_offset;
31         /* objective function has fixed or no integer value */
32         int              is_fixed;
35 #ifdef USE_GMP_FOR_MP
36 #define GBR_type                            mpq_t
37 #define GBR_init(v)                         mpq_init(v)
38 #define GBR_clear(v)                        mpq_clear(v)
39 #define GBR_set(a,b)                        mpq_set(a,b)
40 #define GBR_set_ui(a,b)                     mpq_set_ui(a,b,1)
41 #define GBR_mul(a,b,c)                      mpq_mul(a,b,c)
42 #define GBR_lt(a,b)                         (mpq_cmp(a,b) < 0)
43 #define GBR_is_zero(a)                      (mpq_sgn(a) == 0)
44 #define GBR_numref(a)                       mpq_numref(a)
45 #define GBR_denref(a)                       mpq_denref(a)
46 #define GBR_floor(a,b)                      mpz_fdiv_q(a,GBR_numref(b),GBR_denref(b))
47 #define GBR_ceil(a,b)                       mpz_cdiv_q(a,GBR_numref(b),GBR_denref(b))
48 #endif /* USE_GMP_FOR_MP */
50 #ifdef USE_IMATH_FOR_MP
51 #include <imrat.h>
53 #define GBR_type                            mp_rat
54 #define GBR_init(v)                         v = mp_rat_alloc()
55 #define GBR_clear(v)                        mp_rat_free(v)
56 #define GBR_set(a,b)                        mp_rat_copy(b,a)
57 #define GBR_set_ui(a,b)                     mp_rat_set_uvalue(a,b,1)
58 #define GBR_mul(a,b,c)                      mp_rat_mul(b,c,a)
59 #define GBR_lt(a,b)                         (mp_rat_compare(a,b) < 0)
60 #define GBR_is_zero(a)                      (mp_rat_compare_zero(a) == 0)
61 #define GBR_numref(a)                       mp_rat_numer_ref(a)
62 #define GBR_denref(a)                       mp_rat_denom_ref(a)
63 #define GBR_floor(a,b)                      impz_fdiv_q(a,GBR_numref(b),GBR_denref(b))
64 #define GBR_ceil(a,b)                       impz_cdiv_q(a,GBR_numref(b),GBR_denref(b))
65 #endif /* USE_IMATH_FOR_MP */
67 static struct tab_lp *init_lp(struct isl_tab *tab);
68 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
69 static int solve_lp(struct tab_lp *lp);
70 static void get_obj_val(struct tab_lp* lp, GBR_type *F);
71 static void delete_lp(struct tab_lp *lp);
72 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
73 static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha);
74 static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
75 static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row);
77 #define GBR_LP                              struct tab_lp
78 #define GBR_lp_init(P)                      init_lp(P)
79 #define GBR_lp_set_obj(lp, obj, dim)        set_lp_obj(lp, obj, dim)
80 #define GBR_lp_solve(lp)                    solve_lp(lp)
81 #define GBR_lp_get_obj_val(lp, F)           get_obj_val(lp, F)
82 #define GBR_lp_delete(lp)                   delete_lp(lp)
83 #define GBR_lp_next_row(lp)                 lp->neq
84 #define GBR_lp_add_row(lp, row, dim)        add_lp_row(lp, row, dim)
85 #define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)
86 #define GBR_lp_del_row(lp)                  del_lp_row(lp)
87 #define GBR_lp_is_fixed(lp)                 (lp)->is_fixed
88 #define GBR_lp_cut(lp, obj)                 cut_lp_to_hyperplane(lp, obj)
89 #include "basis_reduction_templ.c"
91 /* Set up a tableau for the Cartesian product of bset with itself.
92  * This could be optimized by first setting up a tableau for bset
93  * and then performing the Cartesian product on the tableau.
94  */
95 static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)
97         unsigned dim;
98         struct isl_tab *prod;
100         if (!tab || !row)
101                 return NULL;
103         dim = tab->n_var;
104         prod = isl_tab_product(tab, tab);
105         if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) {
106                 isl_tab_free(prod);
107                 return NULL;
108         }
109         return prod;
112 static struct tab_lp *init_lp(struct isl_tab *tab)
114         struct tab_lp *lp = NULL;
116         if (!tab)
117                 return NULL;
119         lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);
120         if (!lp)
121                 return NULL;
123         isl_int_init(lp->opt);
124         isl_int_init(lp->opt_denom);
125         isl_int_init(lp->tmp);
126         isl_int_init(lp->tmp2);
128         lp->dim = tab->n_var;
130         lp->ctx = tab->mat->ctx;
131         isl_ctx_ref(lp->ctx);
133         lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
135         lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);
136         if (!lp->row)
137                 goto error;
138         lp->tab = gbr_tab(tab, lp->row);
139         if (!lp->tab)
140                 goto error;
141         lp->con_offset = lp->tab->n_con;
142         lp->obj = NULL;
143         lp->neq = 0;
145         return lp;
146 error:
147         delete_lp(lp);
148         return NULL;
151 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
153         lp->obj = row;
156 static int solve_lp(struct tab_lp *lp)
158         enum isl_lp_result res;
159         unsigned flags = 0;
161         lp->is_fixed = 0;
163         isl_int_set_si(lp->row->el[0], 0);
164         isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);
165         isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);
166         if (lp->neq)
167                 flags = ISL_TAB_SAVE_DUAL;
168         res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
169                           &lp->opt, &lp->opt_denom, flags);
170         isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);
171         if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {
172                 struct isl_vec *sample = isl_tab_get_sample_value(lp->tab);
173                 if (!sample)
174                         return -1;
175                 isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp);
176                 isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2);
177                 isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);
178                 isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);
179                 if (isl_int_ge(lp->tmp, lp->tmp2))
180                         lp->is_fixed = 1;
181                 isl_vec_free(sample);
182         }
183         isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);
184         if (res != isl_lp_ok)
185                 return -1;
186         return 0;
189 /* The current objective function has a fixed (or no) integer value.
190  * Cut the tableau to the hyperplane that fixes this value in
191  * both halves of the tableau.
192  * Return 1 if the resulting tableau is empty.
193  */
194 static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)
196         enum isl_lp_result res;
198         isl_int_set_si(lp->row->el[0], 0);
199         isl_seq_cpy(lp->row->el + 1, row, lp->dim);
200         isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim);
201         res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
202                           &lp->tmp, NULL, 0);
203         if (res != isl_lp_ok)
204                 return -1;
206         isl_int_neg(lp->row->el[0], lp->tmp);
207         if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
208                 return -1;
210         isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim);
211         isl_seq_clr(lp->row->el + 1, lp->dim);
212         if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
213                 return -1;
215         lp->con_offset += 2;
217         return lp->tab->empty;
220 static void get_obj_val(struct tab_lp* lp, GBR_type *F)
222         isl_int_neg(GBR_numref(*F), lp->opt);
223         isl_int_set(GBR_denref(*F), lp->opt_denom);
226 static void delete_lp(struct tab_lp *lp)
228         if (!lp)
229                 return;
231         isl_int_clear(lp->opt);
232         isl_int_clear(lp->opt_denom);
233         isl_int_clear(lp->tmp);
234         isl_int_clear(lp->tmp2);
235         isl_vec_free(lp->row);
236         free(lp->stack);
237         isl_tab_free(lp->tab);
238         isl_ctx_deref(lp->ctx);
239         free(lp);
242 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
244         lp->stack[lp->neq] = isl_tab_snap(lp->tab);
246         isl_int_set_si(lp->row->el[0], 0);
247         isl_seq_cpy(lp->row->el + 1, row, lp->dim);
248         isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);
250         if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0)
251                 return -1;
253         return lp->neq++;
256 static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha)
258         row += lp->con_offset;
259         isl_int_neg(GBR_numref(*alpha), lp->tab->dual->el[1 + row]);
260         isl_int_set(GBR_denref(*alpha), lp->tab->dual->el[0]);
263 static int del_lp_row(struct tab_lp *lp)
265         lp->neq--;
266         return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);