isl_arg_parse: support ISL_ARG_BOOL_ARG flag
[isl.git] / isl_bernstein.c
blob9b8571a509e66fab9ddb6479c5a1521b8fa0829d
1 /*
2 * Copyright 2006-2007 Universiteit Leiden
3 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Copyright 2010 INRIA Saclay
6 * Use of this software is governed by the GNU LGPLv2.1 license
8 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science,
9 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands
10 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A,
11 * B-3001 Leuven, Belgium
12 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
13 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
16 #include <isl/set.h>
17 #include <isl/seq.h>
18 #include <isl_morph.h>
19 #include <isl_factorization.h>
20 #include <isl_vertices_private.h>
21 #include <isl_polynomial_private.h>
22 #include <isl_bernstein.h>
24 struct bernstein_data {
25 enum isl_fold type;
26 isl_qpolynomial *poly;
27 int check_tight;
29 isl_cell *cell;
31 isl_qpolynomial_fold *fold;
32 isl_qpolynomial_fold *fold_tight;
33 isl_pw_qpolynomial_fold *pwf;
34 isl_pw_qpolynomial_fold *pwf_tight;
37 static int vertex_is_integral(__isl_keep isl_basic_set *vertex)
39 unsigned nvar;
40 unsigned nparam;
41 int i;
43 nvar = isl_basic_set_dim(vertex, isl_dim_set);
44 nparam = isl_basic_set_dim(vertex, isl_dim_param);
45 for (i = 0; i < nvar; ++i) {
46 int r = nvar - 1 - i;
47 if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&
48 !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))
49 return 0;
52 return 1;
55 static __isl_give isl_qpolynomial *vertex_coordinate(
56 __isl_keep isl_basic_set *vertex, int i, __isl_take isl_dim *dim)
58 unsigned nvar;
59 unsigned nparam;
60 int r;
61 isl_int denom;
62 isl_qpolynomial *v;
64 nvar = isl_basic_set_dim(vertex, isl_dim_set);
65 nparam = isl_basic_set_dim(vertex, isl_dim_param);
66 r = nvar - 1 - i;
68 isl_int_init(denom);
69 isl_int_set(denom, vertex->eq[r][1 + nparam + i]);
70 isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);
72 if (isl_int_is_pos(denom))
73 isl_seq_neg(vertex->eq[r], vertex->eq[r],
74 1 + isl_basic_set_total_dim(vertex));
75 else
76 isl_int_neg(denom, denom);
78 v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom);
79 isl_int_clear(denom);
81 return v;
82 error:
83 isl_dim_free(dim);
84 isl_int_clear(denom);
85 return NULL;
88 /* Check whether the bound associated to the selection "k" is tight,
89 * which is the case if we select exactly one vertex and if that vertex
90 * is integral for all values of the parameters.
92 static int is_tight(int *k, int n, int d, isl_cell *cell)
94 int i, j;
96 for (i = 0; i < n; ++i) {
97 int v;
98 if (k[i] != d) {
99 if (k[i])
100 return 0;
101 continue;
103 v = cell->ids[n - 1 - i];
104 return vertex_is_integral(cell->vertices->v[v].vertex);
107 return 0;
110 static void add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom,
111 int *k, int n, int d, struct bernstein_data *data)
113 isl_qpolynomial_fold *fold;
115 fold = isl_qpolynomial_fold_alloc(data->type, b);
117 if (data->check_tight && is_tight(k, n, d, data->cell))
118 data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom,
119 data->fold_tight, fold);
120 else
121 data->fold = isl_qpolynomial_fold_fold_on_domain(dom,
122 data->fold, fold);
125 /* Extract the coefficients of the Bernstein base polynomials and store
126 * them in data->fold and data->fold_tight.
128 * In particular, the coefficient of each monomial
129 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding
130 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]!
132 * c[i] contains the coefficient of the selected powers of the first i+1 vars.
133 * multinom[i] contains the partial multinomial coefficient.
135 static void extract_coefficients(isl_qpolynomial *poly,
136 __isl_keep isl_set *dom, struct bernstein_data *data)
138 int i;
139 int d;
140 int n;
141 isl_ctx *ctx;
142 isl_qpolynomial **c = NULL;
143 int *k = NULL;
144 int *left = NULL;
145 isl_vec *multinom = NULL;
147 if (!poly)
148 return;
150 ctx = isl_qpolynomial_get_ctx(poly);
151 n = isl_qpolynomial_dim(poly, isl_dim_set);
152 d = isl_qpolynomial_degree(poly);
153 isl_assert(ctx, n >= 2, return);
155 c = isl_calloc_array(ctx, isl_qpolynomial *, n);
156 k = isl_alloc_array(ctx, int, n);
157 left = isl_alloc_array(ctx, int, n);
158 multinom = isl_vec_alloc(ctx, n);
159 if (!c || !k || !left || !multinom)
160 goto error;
162 isl_int_set_si(multinom->el[0], 1);
163 for (k[0] = d; k[0] >= 0; --k[0]) {
164 int i = 1;
165 isl_qpolynomial_free(c[0]);
166 c[0] = isl_qpolynomial_coeff(poly, isl_dim_set, n - 1, k[0]);
167 left[0] = d - k[0];
168 k[1] = -1;
169 isl_int_set(multinom->el[1], multinom->el[0]);
170 while (i > 0) {
171 if (i == n - 1) {
172 int j;
173 isl_dim *dim;
174 isl_qpolynomial *b;
175 isl_qpolynomial *f;
176 for (j = 2; j <= left[i - 1]; ++j)
177 isl_int_divexact_ui(multinom->el[i],
178 multinom->el[i], j);
179 b = isl_qpolynomial_coeff(c[i - 1], isl_dim_set,
180 n - 1 - i, left[i - 1]);
181 b = isl_qpolynomial_drop_dims(b, isl_dim_set,
182 0, n);
183 dim = isl_qpolynomial_get_dim(b);
184 f = isl_qpolynomial_rat_cst(dim, ctx->one,
185 multinom->el[i]);
186 b = isl_qpolynomial_mul(b, f);
187 k[n - 1] = left[n - 2];
188 add_fold(b, dom, k, n, d, data);
189 --i;
190 continue;
192 if (k[i] >= left[i - 1]) {
193 --i;
194 continue;
196 ++k[i];
197 if (k[i])
198 isl_int_divexact_ui(multinom->el[i],
199 multinom->el[i], k[i]);
200 isl_qpolynomial_free(c[i]);
201 c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_set,
202 n - 1 - i, k[i]);
203 left[i] = left[i - 1] - k[i];
204 k[i + 1] = -1;
205 isl_int_set(multinom->el[i + 1], multinom->el[i]);
206 ++i;
208 isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);
211 for (i = 0; i < n; ++i)
212 isl_qpolynomial_free(c[i]);
214 isl_vec_free(multinom);
215 free(left);
216 free(k);
217 free(c);
218 return;
219 error:
220 isl_vec_free(multinom);
221 free(left);
222 free(k);
223 if (c)
224 for (i = 0; i < n; ++i)
225 isl_qpolynomial_free(c[i]);
226 free(c);
227 return;
230 /* Perform bernstein expansion on the parametric vertices that are active
231 * on "cell".
233 * data->poly has been homogenized in the calling function.
235 * We plug in the barycentric coordinates for the set variables
237 * \vec x = \sum_i \alpha_i v_i(\vec p)
239 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.
240 * Next, we extract the coefficients of the Bernstein base polynomials.
242 static int bernstein_coefficients_cell(__isl_take isl_cell *cell, void *user)
244 int i, j;
245 struct bernstein_data *data = (struct bernstein_data *)user;
246 isl_dim *dim_param;
247 isl_dim *dim_dst;
248 isl_qpolynomial *poly = data->poly;
249 unsigned nvar;
250 int n_vertices;
251 isl_qpolynomial **subs;
252 isl_pw_qpolynomial_fold *pwf;
253 isl_set *dom;
254 isl_ctx *ctx;
256 nvar = isl_qpolynomial_dim(poly, isl_dim_set) - 1;
257 n_vertices = cell->n_vertices;
259 ctx = isl_qpolynomial_get_ctx(poly);
260 if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate)
261 return isl_cell_foreach_simplex(cell,
262 &bernstein_coefficients_cell, user);
264 subs = isl_alloc_array(data->poly->dim->ctx, isl_qpolynomial *,
265 1 + nvar);
266 if (!subs)
267 goto error;
269 dim_param = isl_basic_set_get_dim(cell->dom);
270 dim_dst = isl_qpolynomial_get_dim(poly);
271 dim_dst = isl_dim_add(dim_dst, isl_dim_set, n_vertices);
273 for (i = 0; i < 1 + nvar; ++i)
274 subs[i] = isl_qpolynomial_zero(isl_dim_copy(dim_dst));
276 for (i = 0; i < n_vertices; ++i) {
277 isl_qpolynomial *c;
278 c = isl_qpolynomial_var(isl_dim_copy(dim_dst), isl_dim_set,
279 1 + nvar + i);
280 for (j = 0; j < nvar; ++j) {
281 int k = cell->ids[i];
282 isl_qpolynomial *v;
283 v = vertex_coordinate(cell->vertices->v[k].vertex, j,
284 isl_dim_copy(dim_param));
285 v = isl_qpolynomial_add_dims(v, isl_dim_set,
286 1 + nvar + n_vertices);
287 v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c));
288 subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v);
290 subs[0] = isl_qpolynomial_add(subs[0], c);
292 isl_dim_free(dim_dst);
294 poly = isl_qpolynomial_copy(poly);
296 poly = isl_qpolynomial_add_dims(poly, isl_dim_set, n_vertices);
297 poly = isl_qpolynomial_substitute(poly, isl_dim_set, 0, 1 + nvar, subs);
298 poly = isl_qpolynomial_drop_dims(poly, isl_dim_set, 0, 1 + nvar);
300 data->cell = cell;
301 dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom));
302 data->fold = isl_qpolynomial_fold_empty(data->type, isl_dim_copy(dim_param));
303 data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param);
304 extract_coefficients(poly, dom, data);
306 pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),
307 data->fold);
308 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);
309 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight);
310 data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf);
312 isl_qpolynomial_free(poly);
313 isl_cell_free(cell);
314 for (i = 0; i < 1 + nvar; ++i)
315 isl_qpolynomial_free(subs[i]);
316 free(subs);
317 return 0;
318 error:
319 isl_cell_free(cell);
320 return -1;
323 /* Base case of applying bernstein expansion.
325 * We compute the chamber decomposition of the parametric polytope "bset"
326 * and then perform bernstein expansion on the parametric vertices
327 * that are active on each chamber.
329 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base(
330 __isl_take isl_basic_set *bset,
331 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
333 unsigned nvar;
334 isl_dim *dim;
335 isl_pw_qpolynomial_fold *pwf;
336 isl_vertices *vertices;
337 int covers;
339 nvar = isl_basic_set_dim(bset, isl_dim_set);
340 if (nvar == 0) {
341 isl_set *dom;
342 isl_qpolynomial_fold *fold;
343 fold = isl_qpolynomial_fold_alloc(data->type, poly);
344 dom = isl_set_from_basic_set(bset);
345 if (tight)
346 *tight = 1;
347 return isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
350 if (isl_qpolynomial_is_zero(poly)) {
351 isl_set *dom;
352 isl_qpolynomial_fold *fold;
353 fold = isl_qpolynomial_fold_alloc(data->type, poly);
354 dom = isl_set_from_basic_set(bset);
355 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
356 if (tight)
357 *tight = 1;
358 return isl_pw_qpolynomial_fold_drop_dims(pwf,
359 isl_dim_set, 0, nvar);
362 dim = isl_basic_set_get_dim(bset);
363 dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
364 data->pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim), data->type);
365 data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type);
366 data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly));
367 vertices = isl_basic_set_compute_vertices(bset);
368 isl_vertices_foreach_disjoint_cell(vertices,
369 &bernstein_coefficients_cell, data);
370 isl_vertices_free(vertices);
371 isl_qpolynomial_free(data->poly);
373 isl_basic_set_free(bset);
374 isl_qpolynomial_free(poly);
376 covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);
377 if (covers < 0)
378 goto error;
380 if (tight)
381 *tight = covers;
383 if (covers) {
384 isl_pw_qpolynomial_fold_free(data->pwf);
385 return data->pwf_tight;
388 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);
390 return data->pwf;
391 error:
392 isl_pw_qpolynomial_fold_free(data->pwf_tight);
393 isl_pw_qpolynomial_fold_free(data->pwf);
394 return NULL;
397 /* Apply bernstein expansion recursively by working in on len[i]
398 * set variables at a time, with i ranging from n_group - 1 to 0.
400 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive(
401 __isl_take isl_pw_qpolynomial *pwqp,
402 int n_group, int *len, struct bernstein_data *data, int *tight)
404 int i;
405 unsigned nparam;
406 unsigned nvar;
407 isl_pw_qpolynomial_fold *pwf;
409 if (!pwqp)
410 return NULL;
412 nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
413 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_set);
415 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,
416 isl_dim_set, 0, nvar - len[n_group - 1]);
417 pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight);
419 for (i = n_group - 2; i >= 0; --i) {
420 nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);
421 pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_set, 0,
422 isl_dim_param, nparam - len[i], len[i]);
423 if (tight && !*tight)
424 tight = NULL;
425 pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);
428 return pwf;
431 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors(
432 __isl_take isl_basic_set *bset,
433 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
435 isl_factorizer *f;
436 isl_set *set;
437 isl_pw_qpolynomial *pwqp;
438 isl_pw_qpolynomial_fold *pwf;
440 f = isl_basic_set_factorizer(bset);
441 if (!f)
442 goto error;
443 if (f->n_group == 0) {
444 isl_factorizer_free(f);
445 return bernstein_coefficients_base(bset, poly, data, tight);
448 set = isl_set_from_basic_set(bset);
449 pwqp = isl_pw_qpolynomial_alloc(set, poly);
450 pwqp = isl_pw_qpolynomial_morph(pwqp, isl_morph_copy(f->morph));
452 pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,
453 tight);
455 isl_factorizer_free(f);
457 return pwf;
458 error:
459 isl_basic_set_free(bset);
460 isl_qpolynomial_free(poly);
461 return NULL;
464 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive(
465 __isl_take isl_basic_set *bset,
466 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
468 int i;
469 int *len;
470 unsigned nvar;
471 isl_pw_qpolynomial_fold *pwf;
472 isl_set *set;
473 isl_pw_qpolynomial *pwqp;
475 if (!bset || !poly)
476 goto error;
478 nvar = isl_basic_set_dim(bset, isl_dim_set);
480 len = isl_alloc_array(bset->ctx, int, nvar);
481 if (!len)
482 goto error;
484 for (i = 0; i < nvar; ++i)
485 len[i] = 1;
487 set = isl_set_from_basic_set(bset);
488 pwqp = isl_pw_qpolynomial_alloc(set, poly);
490 pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);
492 free(len);
494 return pwf;
495 error:
496 isl_basic_set_free(bset);
497 isl_qpolynomial_free(poly);
498 return NULL;
501 /* Compute a bound on the polynomial defined over the parametric polytope
502 * using bernstein expansion and store the result
503 * in bound->pwf and bound->pwf_tight.
505 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if
506 * the polytope can be factorized and apply bernstein expansion recursively
507 * on the factors.
508 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply
509 * bernstein expansion recursively on each dimension.
510 * Otherwise, we apply bernstein expansion on the entire polytope.
512 int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set *bset,
513 __isl_take isl_qpolynomial *poly, struct isl_bound *bound)
515 struct bernstein_data data;
516 isl_pw_qpolynomial_fold *pwf;
517 unsigned nvar;
518 int tight = 0;
519 int *tp = bound->check_tight ? &tight : NULL;
521 if (!bset || !poly)
522 goto error;
524 data.type = bound->type;
525 data.check_tight = bound->check_tight;
527 nvar = isl_basic_set_dim(bset, isl_dim_set);
529 if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)
530 pwf = bernstein_coefficients_factors(bset, poly, &data, tp);
531 else if (nvar > 1 &&
532 (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))
533 pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);
534 else
535 pwf = bernstein_coefficients_base(bset, poly, &data, tp);
537 if (tight)
538 bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf);
539 else
540 bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf);
542 return 0;
543 error:
544 isl_basic_set_free(bset);
545 isl_qpolynomial_free(poly);
546 return -1;