isl_space_from_range: rename "dim" argument to "space"
[isl.git] / isl_bernstein.c
blob2bc2974764fdb8aa1a449bde9e02d7c4c59eb4b6
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 MIT 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_ctx_private.h>
17 #include <isl_map_private.h>
18 #include <isl/set.h>
19 #include <isl_seq.h>
20 #include <isl_morph.h>
21 #include <isl_factorization.h>
22 #include <isl_vertices_private.h>
23 #include <isl_polynomial_private.h>
24 #include <isl_options_private.h>
25 #include <isl_vec_private.h>
26 #include <isl_bernstein.h>
28 struct bernstein_data {
29 enum isl_fold type;
30 isl_qpolynomial *poly;
31 int check_tight;
33 isl_cell *cell;
35 isl_qpolynomial_fold *fold;
36 isl_qpolynomial_fold *fold_tight;
37 isl_pw_qpolynomial_fold *pwf;
38 isl_pw_qpolynomial_fold *pwf_tight;
41 static int vertex_is_integral(__isl_keep isl_basic_set *vertex)
43 unsigned nvar;
44 unsigned nparam;
45 int i;
47 nvar = isl_basic_set_dim(vertex, isl_dim_set);
48 nparam = isl_basic_set_dim(vertex, isl_dim_param);
49 for (i = 0; i < nvar; ++i) {
50 int r = nvar - 1 - i;
51 if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&
52 !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))
53 return 0;
56 return 1;
59 static __isl_give isl_qpolynomial *vertex_coordinate(
60 __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *space)
62 unsigned nvar;
63 unsigned nparam;
64 int r;
65 isl_int denom;
66 isl_qpolynomial *v;
68 nvar = isl_basic_set_dim(vertex, isl_dim_set);
69 nparam = isl_basic_set_dim(vertex, isl_dim_param);
70 r = nvar - 1 - i;
72 isl_int_init(denom);
73 isl_int_set(denom, vertex->eq[r][1 + nparam + i]);
74 isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);
76 if (isl_int_is_pos(denom))
77 isl_seq_neg(vertex->eq[r], vertex->eq[r],
78 1 + isl_basic_set_total_dim(vertex));
79 else
80 isl_int_neg(denom, denom);
82 v = isl_qpolynomial_from_affine(space, vertex->eq[r], denom);
83 isl_int_clear(denom);
85 return v;
86 error:
87 isl_space_free(space);
88 isl_int_clear(denom);
89 return NULL;
92 /* Check whether the bound associated to the selection "k" is tight,
93 * which is the case if we select exactly one vertex and if that vertex
94 * is integral for all values of the parameters.
96 static int is_tight(int *k, int n, int d, isl_cell *cell)
98 int i;
100 for (i = 0; i < n; ++i) {
101 int v;
102 if (k[i] != d) {
103 if (k[i])
104 return 0;
105 continue;
107 v = cell->ids[n - 1 - i];
108 return vertex_is_integral(cell->vertices->v[v].vertex);
111 return 0;
114 static void add_fold(__isl_take isl_qpolynomial *b, __isl_keep isl_set *dom,
115 int *k, int n, int d, struct bernstein_data *data)
117 isl_qpolynomial_fold *fold;
119 fold = isl_qpolynomial_fold_alloc(data->type, b);
121 if (data->check_tight && is_tight(k, n, d, data->cell))
122 data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom,
123 data->fold_tight, fold);
124 else
125 data->fold = isl_qpolynomial_fold_fold_on_domain(dom,
126 data->fold, fold);
129 /* Extract the coefficients of the Bernstein base polynomials and store
130 * them in data->fold and data->fold_tight.
132 * In particular, the coefficient of each monomial
133 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding
134 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]!
136 * c[i] contains the coefficient of the selected powers of the first i+1 vars.
137 * multinom[i] contains the partial multinomial coefficient.
139 static isl_stat extract_coefficients(isl_qpolynomial *poly,
140 __isl_keep isl_set *dom, struct bernstein_data *data)
142 int i;
143 int d;
144 int n;
145 isl_ctx *ctx;
146 isl_qpolynomial **c = NULL;
147 int *k = NULL;
148 int *left = NULL;
149 isl_vec *multinom = NULL;
151 if (!poly)
152 return isl_stat_error;
154 ctx = isl_qpolynomial_get_ctx(poly);
155 n = isl_qpolynomial_dim(poly, isl_dim_in);
156 d = isl_qpolynomial_degree(poly);
157 isl_assert(ctx, n >= 2, return isl_stat_error);
159 c = isl_calloc_array(ctx, isl_qpolynomial *, n);
160 k = isl_alloc_array(ctx, int, n);
161 left = isl_alloc_array(ctx, int, n);
162 multinom = isl_vec_alloc(ctx, n);
163 if (!c || !k || !left || !multinom)
164 goto error;
166 isl_int_set_si(multinom->el[0], 1);
167 for (k[0] = d; k[0] >= 0; --k[0]) {
168 int i = 1;
169 isl_qpolynomial_free(c[0]);
170 c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]);
171 left[0] = d - k[0];
172 k[1] = -1;
173 isl_int_set(multinom->el[1], multinom->el[0]);
174 while (i > 0) {
175 if (i == n - 1) {
176 int j;
177 isl_space *dim;
178 isl_qpolynomial *b;
179 isl_qpolynomial *f;
180 for (j = 2; j <= left[i - 1]; ++j)
181 isl_int_divexact_ui(multinom->el[i],
182 multinom->el[i], j);
183 b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,
184 n - 1 - i, left[i - 1]);
185 b = isl_qpolynomial_project_domain_on_params(b);
186 dim = isl_qpolynomial_get_domain_space(b);
187 f = isl_qpolynomial_rat_cst_on_domain(dim, ctx->one,
188 multinom->el[i]);
189 b = isl_qpolynomial_mul(b, f);
190 k[n - 1] = left[n - 2];
191 add_fold(b, dom, k, n, d, data);
192 --i;
193 continue;
195 if (k[i] >= left[i - 1]) {
196 --i;
197 continue;
199 ++k[i];
200 if (k[i])
201 isl_int_divexact_ui(multinom->el[i],
202 multinom->el[i], k[i]);
203 isl_qpolynomial_free(c[i]);
204 c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,
205 n - 1 - i, k[i]);
206 left[i] = left[i - 1] - k[i];
207 k[i + 1] = -1;
208 isl_int_set(multinom->el[i + 1], multinom->el[i]);
209 ++i;
211 isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);
214 for (i = 0; i < n; ++i)
215 isl_qpolynomial_free(c[i]);
217 isl_vec_free(multinom);
218 free(left);
219 free(k);
220 free(c);
221 return isl_stat_ok;
222 error:
223 isl_vec_free(multinom);
224 free(left);
225 free(k);
226 if (c)
227 for (i = 0; i < n; ++i)
228 isl_qpolynomial_free(c[i]);
229 free(c);
230 return isl_stat_error;
233 /* Perform bernstein expansion on the parametric vertices that are active
234 * on "cell".
236 * data->poly has been homogenized in the calling function.
238 * We plug in the barycentric coordinates for the set variables
240 * \vec x = \sum_i \alpha_i v_i(\vec p)
242 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.
243 * Next, we extract the coefficients of the Bernstein base polynomials.
245 static isl_stat bernstein_coefficients_cell(__isl_take isl_cell *cell,
246 void *user)
248 int i, j;
249 struct bernstein_data *data = (struct bernstein_data *)user;
250 isl_space *dim_param;
251 isl_space *dim_dst;
252 isl_qpolynomial *poly = data->poly;
253 unsigned nvar;
254 int n_vertices;
255 isl_qpolynomial **subs;
256 isl_pw_qpolynomial_fold *pwf;
257 isl_set *dom;
258 isl_ctx *ctx;
260 if (!poly)
261 goto error;
263 nvar = isl_qpolynomial_dim(poly, isl_dim_in) - 1;
264 n_vertices = cell->n_vertices;
266 ctx = isl_qpolynomial_get_ctx(poly);
267 if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate)
268 return isl_cell_foreach_simplex(cell,
269 &bernstein_coefficients_cell, user);
271 subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar);
272 if (!subs)
273 goto error;
275 dim_param = isl_basic_set_get_space(cell->dom);
276 dim_dst = isl_qpolynomial_get_domain_space(poly);
277 dim_dst = isl_space_add_dims(dim_dst, isl_dim_set, n_vertices);
279 for (i = 0; i < 1 + nvar; ++i)
280 subs[i] = isl_qpolynomial_zero_on_domain(isl_space_copy(dim_dst));
282 for (i = 0; i < n_vertices; ++i) {
283 isl_qpolynomial *c;
284 c = isl_qpolynomial_var_on_domain(isl_space_copy(dim_dst), isl_dim_set,
285 1 + nvar + i);
286 for (j = 0; j < nvar; ++j) {
287 int k = cell->ids[i];
288 isl_qpolynomial *v;
289 v = vertex_coordinate(cell->vertices->v[k].vertex, j,
290 isl_space_copy(dim_param));
291 v = isl_qpolynomial_add_dims(v, isl_dim_in,
292 1 + nvar + n_vertices);
293 v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c));
294 subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v);
296 subs[0] = isl_qpolynomial_add(subs[0], c);
298 isl_space_free(dim_dst);
300 poly = isl_qpolynomial_copy(poly);
302 poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices);
303 poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs);
304 poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar);
306 data->cell = cell;
307 dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom));
308 data->fold = isl_qpolynomial_fold_empty(data->type, isl_space_copy(dim_param));
309 data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param);
310 if (extract_coefficients(poly, dom, data) < 0) {
311 data->fold = isl_qpolynomial_fold_free(data->fold);
312 data->fold_tight = isl_qpolynomial_fold_free(data->fold_tight);
315 pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),
316 data->fold);
317 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);
318 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight);
319 data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf);
321 isl_qpolynomial_free(poly);
322 isl_cell_free(cell);
323 for (i = 0; i < 1 + nvar; ++i)
324 isl_qpolynomial_free(subs[i]);
325 free(subs);
326 return isl_stat_ok;
327 error:
328 isl_cell_free(cell);
329 return isl_stat_error;
332 /* Base case of applying bernstein expansion.
334 * We compute the chamber decomposition of the parametric polytope "bset"
335 * and then perform bernstein expansion on the parametric vertices
336 * that are active on each chamber.
338 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base(
339 __isl_take isl_basic_set *bset,
340 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
342 unsigned nvar;
343 isl_space *dim;
344 isl_pw_qpolynomial_fold *pwf;
345 isl_vertices *vertices;
346 int covers;
348 nvar = isl_basic_set_dim(bset, isl_dim_set);
349 if (nvar == 0) {
350 isl_set *dom;
351 isl_qpolynomial_fold *fold;
353 fold = isl_qpolynomial_fold_alloc(data->type, poly);
354 dom = isl_set_from_basic_set(bset);
355 if (tight)
356 *tight = 1;
357 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
358 return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);
361 if (isl_qpolynomial_is_zero(poly)) {
362 isl_set *dom;
363 isl_qpolynomial_fold *fold;
364 fold = isl_qpolynomial_fold_alloc(data->type, poly);
365 dom = isl_set_from_basic_set(bset);
366 pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
367 if (tight)
368 *tight = 1;
369 return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);
372 dim = isl_basic_set_get_space(bset);
373 dim = isl_space_params(dim);
374 dim = isl_space_from_domain(dim);
375 dim = isl_space_add_dims(dim, isl_dim_set, 1);
376 data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(dim), data->type);
377 data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type);
378 data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly));
379 vertices = isl_basic_set_compute_vertices(bset);
380 if (isl_vertices_foreach_disjoint_cell(vertices,
381 &bernstein_coefficients_cell, data) < 0)
382 data->pwf = isl_pw_qpolynomial_fold_free(data->pwf);
383 isl_vertices_free(vertices);
384 isl_qpolynomial_free(data->poly);
386 isl_basic_set_free(bset);
387 isl_qpolynomial_free(poly);
389 covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);
390 if (covers < 0)
391 goto error;
393 if (tight)
394 *tight = covers;
396 if (covers) {
397 isl_pw_qpolynomial_fold_free(data->pwf);
398 return data->pwf_tight;
401 data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);
403 return data->pwf;
404 error:
405 isl_pw_qpolynomial_fold_free(data->pwf_tight);
406 isl_pw_qpolynomial_fold_free(data->pwf);
407 return NULL;
410 /* Apply bernstein expansion recursively by working in on len[i]
411 * set variables at a time, with i ranging from n_group - 1 to 0.
413 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive(
414 __isl_take isl_pw_qpolynomial *pwqp,
415 int n_group, int *len, struct bernstein_data *data, int *tight)
417 int i;
418 unsigned nparam;
419 unsigned nvar;
420 isl_pw_qpolynomial_fold *pwf;
422 if (!pwqp)
423 return NULL;
425 nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
426 nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in);
428 pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,
429 isl_dim_in, 0, nvar - len[n_group - 1]);
430 pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight);
432 for (i = n_group - 2; i >= 0; --i) {
433 nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);
434 pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0,
435 isl_dim_param, nparam - len[i], len[i]);
436 if (tight && !*tight)
437 tight = NULL;
438 pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);
441 return pwf;
444 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors(
445 __isl_take isl_basic_set *bset,
446 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
448 isl_factorizer *f;
449 isl_set *set;
450 isl_pw_qpolynomial *pwqp;
451 isl_pw_qpolynomial_fold *pwf;
453 f = isl_basic_set_factorizer(bset);
454 if (!f)
455 goto error;
456 if (f->n_group == 0) {
457 isl_factorizer_free(f);
458 return bernstein_coefficients_base(bset, poly, data, tight);
461 set = isl_set_from_basic_set(bset);
462 pwqp = isl_pw_qpolynomial_alloc(set, poly);
463 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph));
465 pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,
466 tight);
468 isl_factorizer_free(f);
470 return pwf;
471 error:
472 isl_basic_set_free(bset);
473 isl_qpolynomial_free(poly);
474 return NULL;
477 static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive(
478 __isl_take isl_basic_set *bset,
479 __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
481 int i;
482 int *len;
483 unsigned nvar;
484 isl_pw_qpolynomial_fold *pwf;
485 isl_set *set;
486 isl_pw_qpolynomial *pwqp;
488 if (!bset || !poly)
489 goto error;
491 nvar = isl_basic_set_dim(bset, isl_dim_set);
493 len = isl_alloc_array(bset->ctx, int, nvar);
494 if (nvar && !len)
495 goto error;
497 for (i = 0; i < nvar; ++i)
498 len[i] = 1;
500 set = isl_set_from_basic_set(bset);
501 pwqp = isl_pw_qpolynomial_alloc(set, poly);
503 pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);
505 free(len);
507 return pwf;
508 error:
509 isl_basic_set_free(bset);
510 isl_qpolynomial_free(poly);
511 return NULL;
514 /* Compute a bound on the polynomial defined over the parametric polytope
515 * using bernstein expansion and store the result
516 * in bound->pwf and bound->pwf_tight.
518 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if
519 * the polytope can be factorized and apply bernstein expansion recursively
520 * on the factors.
521 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply
522 * bernstein expansion recursively on each dimension.
523 * Otherwise, we apply bernstein expansion on the entire polytope.
525 isl_stat isl_qpolynomial_bound_on_domain_bernstein(
526 __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *poly,
527 struct isl_bound *bound)
529 struct bernstein_data data;
530 isl_pw_qpolynomial_fold *pwf;
531 unsigned nvar;
532 int tight = 0;
533 int *tp = bound->check_tight ? &tight : NULL;
535 if (!bset || !poly)
536 goto error;
538 data.type = bound->type;
539 data.check_tight = bound->check_tight;
541 nvar = isl_basic_set_dim(bset, isl_dim_set);
543 if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)
544 pwf = bernstein_coefficients_factors(bset, poly, &data, tp);
545 else if (nvar > 1 &&
546 (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))
547 pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);
548 else
549 pwf = bernstein_coefficients_base(bset, poly, &data, tp);
551 if (tight)
552 bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf);
553 else
554 bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf);
556 return isl_stat_ok;
557 error:
558 isl_basic_set_free(bset);
559 isl_qpolynomial_free(poly);
560 return isl_stat_error;