Line data Source code
1 : /*
2 : * Copyright 2006-2007 Universiteit Leiden
3 : * Copyright 2008-2009 Katholieke Universiteit Leuven
4 : * Copyright 2010 INRIA Saclay
5 : *
6 : * Use of this software is governed by the MIT license
7 : *
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
14 : */
15 :
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>
27 :
28 : struct bernstein_data {
29 : enum isl_fold type;
30 : isl_qpolynomial *poly;
31 : int check_tight;
32 :
33 : isl_cell *cell;
34 :
35 : isl_qpolynomial_fold *fold;
36 : isl_qpolynomial_fold *fold_tight;
37 : isl_pw_qpolynomial_fold *pwf;
38 : isl_pw_qpolynomial_fold *pwf_tight;
39 : };
40 :
41 0 : static int vertex_is_integral(__isl_keep isl_basic_set *vertex)
42 : {
43 : unsigned nvar;
44 : unsigned nparam;
45 : int i;
46 :
47 0 : nvar = isl_basic_set_dim(vertex, isl_dim_set);
48 0 : nparam = isl_basic_set_dim(vertex, isl_dim_param);
49 0 : for (i = 0; i < nvar; ++i) {
50 0 : int r = nvar - 1 - i;
51 0 : if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&
52 0 : !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))
53 0 : return 0;
54 : }
55 :
56 0 : return 1;
57 : }
58 :
59 0 : static __isl_give isl_qpolynomial *vertex_coordinate(
60 : __isl_keep isl_basic_set *vertex, int i, __isl_take isl_space *dim)
61 : {
62 : unsigned nvar;
63 : unsigned nparam;
64 : int r;
65 : isl_int denom;
66 : isl_qpolynomial *v;
67 :
68 0 : nvar = isl_basic_set_dim(vertex, isl_dim_set);
69 0 : nparam = isl_basic_set_dim(vertex, isl_dim_param);
70 0 : r = nvar - 1 - i;
71 :
72 0 : isl_int_init(denom);
73 0 : isl_int_set(denom, vertex->eq[r][1 + nparam + i]);
74 0 : isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);
75 :
76 0 : if (isl_int_is_pos(denom))
77 0 : isl_seq_neg(vertex->eq[r], vertex->eq[r],
78 0 : 1 + isl_basic_set_total_dim(vertex));
79 : else
80 0 : isl_int_neg(denom, denom);
81 :
82 0 : v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom);
83 0 : isl_int_clear(denom);
84 :
85 0 : return v;
86 : error:
87 0 : isl_space_free(dim);
88 0 : isl_int_clear(denom);
89 0 : return NULL;
90 : }
91 :
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.
95 : */
96 0 : static int is_tight(int *k, int n, int d, isl_cell *cell)
97 : {
98 : int i;
99 :
100 0 : for (i = 0; i < n; ++i) {
101 : int v;
102 0 : if (k[i] != d) {
103 0 : if (k[i])
104 0 : return 0;
105 0 : continue;
106 : }
107 0 : v = cell->ids[n - 1 - i];
108 0 : return vertex_is_integral(cell->vertices->v[v].vertex);
109 : }
110 :
111 0 : return 0;
112 : }
113 :
114 0 : 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)
116 : {
117 : isl_qpolynomial_fold *fold;
118 :
119 0 : fold = isl_qpolynomial_fold_alloc(data->type, b);
120 :
121 0 : if (data->check_tight && is_tight(k, n, d, data->cell))
122 0 : data->fold_tight = isl_qpolynomial_fold_fold_on_domain(dom,
123 : data->fold_tight, fold);
124 : else
125 0 : data->fold = isl_qpolynomial_fold_fold_on_domain(dom,
126 : data->fold, fold);
127 0 : }
128 :
129 : /* Extract the coefficients of the Bernstein base polynomials and store
130 : * them in data->fold and data->fold_tight.
131 : *
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]!
135 : *
136 : * c[i] contains the coefficient of the selected powers of the first i+1 vars.
137 : * multinom[i] contains the partial multinomial coefficient.
138 : */
139 0 : static void extract_coefficients(isl_qpolynomial *poly,
140 : __isl_keep isl_set *dom, struct bernstein_data *data)
141 : {
142 : int i;
143 : int d;
144 : int n;
145 : isl_ctx *ctx;
146 0 : isl_qpolynomial **c = NULL;
147 0 : int *k = NULL;
148 0 : int *left = NULL;
149 0 : isl_vec *multinom = NULL;
150 :
151 0 : if (!poly)
152 0 : return;
153 :
154 0 : ctx = isl_qpolynomial_get_ctx(poly);
155 0 : n = isl_qpolynomial_dim(poly, isl_dim_in);
156 0 : d = isl_qpolynomial_degree(poly);
157 0 : isl_assert(ctx, n >= 2, return);
158 :
159 0 : c = isl_calloc_array(ctx, isl_qpolynomial *, n);
160 0 : k = isl_alloc_array(ctx, int, n);
161 0 : left = isl_alloc_array(ctx, int, n);
162 0 : multinom = isl_vec_alloc(ctx, n);
163 0 : if (!c || !k || !left || !multinom)
164 : goto error;
165 :
166 0 : isl_int_set_si(multinom->el[0], 1);
167 0 : for (k[0] = d; k[0] >= 0; --k[0]) {
168 0 : int i = 1;
169 0 : isl_qpolynomial_free(c[0]);
170 0 : c[0] = isl_qpolynomial_coeff(poly, isl_dim_in, n - 1, k[0]);
171 0 : left[0] = d - k[0];
172 0 : k[1] = -1;
173 0 : isl_int_set(multinom->el[1], multinom->el[0]);
174 0 : while (i > 0) {
175 0 : if (i == n - 1) {
176 : int j;
177 : isl_space *dim;
178 : isl_qpolynomial *b;
179 : isl_qpolynomial *f;
180 0 : for (j = 2; j <= left[i - 1]; ++j)
181 0 : isl_int_divexact_ui(multinom->el[i],
182 : multinom->el[i], j);
183 0 : b = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,
184 0 : n - 1 - i, left[i - 1]);
185 0 : b = isl_qpolynomial_project_domain_on_params(b);
186 0 : dim = isl_qpolynomial_get_domain_space(b);
187 0 : f = isl_qpolynomial_rat_cst_on_domain(dim, ctx->one,
188 0 : multinom->el[i]);
189 0 : b = isl_qpolynomial_mul(b, f);
190 0 : k[n - 1] = left[n - 2];
191 0 : add_fold(b, dom, k, n, d, data);
192 0 : --i;
193 0 : continue;
194 : }
195 0 : if (k[i] >= left[i - 1]) {
196 0 : --i;
197 0 : continue;
198 : }
199 0 : ++k[i];
200 0 : if (k[i])
201 0 : isl_int_divexact_ui(multinom->el[i],
202 : multinom->el[i], k[i]);
203 0 : isl_qpolynomial_free(c[i]);
204 0 : c[i] = isl_qpolynomial_coeff(c[i - 1], isl_dim_in,
205 0 : n - 1 - i, k[i]);
206 0 : left[i] = left[i - 1] - k[i];
207 0 : k[i + 1] = -1;
208 0 : isl_int_set(multinom->el[i + 1], multinom->el[i]);
209 0 : ++i;
210 : }
211 0 : isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);
212 : }
213 :
214 0 : for (i = 0; i < n; ++i)
215 0 : isl_qpolynomial_free(c[i]);
216 :
217 0 : isl_vec_free(multinom);
218 0 : free(left);
219 0 : free(k);
220 0 : free(c);
221 0 : return;
222 : error:
223 0 : isl_vec_free(multinom);
224 0 : free(left);
225 0 : free(k);
226 0 : if (c)
227 0 : for (i = 0; i < n; ++i)
228 0 : isl_qpolynomial_free(c[i]);
229 0 : free(c);
230 0 : return;
231 : }
232 :
233 : /* Perform bernstein expansion on the parametric vertices that are active
234 : * on "cell".
235 : *
236 : * data->poly has been homogenized in the calling function.
237 : *
238 : * We plug in the barycentric coordinates for the set variables
239 : *
240 : * \vec x = \sum_i \alpha_i v_i(\vec p)
241 : *
242 : * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.
243 : * Next, we extract the coefficients of the Bernstein base polynomials.
244 : */
245 0 : static isl_stat bernstein_coefficients_cell(__isl_take isl_cell *cell,
246 : void *user)
247 : {
248 : int i, j;
249 0 : struct bernstein_data *data = (struct bernstein_data *)user;
250 : isl_space *dim_param;
251 : isl_space *dim_dst;
252 0 : 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;
259 :
260 0 : if (!poly)
261 0 : goto error;
262 :
263 0 : nvar = isl_qpolynomial_dim(poly, isl_dim_in) - 1;
264 0 : n_vertices = cell->n_vertices;
265 :
266 0 : ctx = isl_qpolynomial_get_ctx(poly);
267 0 : if (n_vertices > nvar + 1 && ctx->opt->bernstein_triangulate)
268 0 : return isl_cell_foreach_simplex(cell,
269 : &bernstein_coefficients_cell, user);
270 :
271 0 : subs = isl_alloc_array(ctx, isl_qpolynomial *, 1 + nvar);
272 0 : if (!subs)
273 0 : goto error;
274 :
275 0 : dim_param = isl_basic_set_get_space(cell->dom);
276 0 : dim_dst = isl_qpolynomial_get_domain_space(poly);
277 0 : dim_dst = isl_space_add_dims(dim_dst, isl_dim_set, n_vertices);
278 :
279 0 : for (i = 0; i < 1 + nvar; ++i)
280 0 : subs[i] = isl_qpolynomial_zero_on_domain(isl_space_copy(dim_dst));
281 :
282 0 : for (i = 0; i < n_vertices; ++i) {
283 : isl_qpolynomial *c;
284 0 : c = isl_qpolynomial_var_on_domain(isl_space_copy(dim_dst), isl_dim_set,
285 0 : 1 + nvar + i);
286 0 : for (j = 0; j < nvar; ++j) {
287 0 : int k = cell->ids[i];
288 : isl_qpolynomial *v;
289 0 : v = vertex_coordinate(cell->vertices->v[k].vertex, j,
290 : isl_space_copy(dim_param));
291 0 : v = isl_qpolynomial_add_dims(v, isl_dim_in,
292 0 : 1 + nvar + n_vertices);
293 0 : v = isl_qpolynomial_mul(v, isl_qpolynomial_copy(c));
294 0 : subs[1 + j] = isl_qpolynomial_add(subs[1 + j], v);
295 : }
296 0 : subs[0] = isl_qpolynomial_add(subs[0], c);
297 : }
298 0 : isl_space_free(dim_dst);
299 :
300 0 : poly = isl_qpolynomial_copy(poly);
301 :
302 0 : poly = isl_qpolynomial_add_dims(poly, isl_dim_in, n_vertices);
303 0 : poly = isl_qpolynomial_substitute(poly, isl_dim_in, 0, 1 + nvar, subs);
304 0 : poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, 1 + nvar);
305 :
306 0 : data->cell = cell;
307 0 : dom = isl_set_from_basic_set(isl_basic_set_copy(cell->dom));
308 0 : data->fold = isl_qpolynomial_fold_empty(data->type, isl_space_copy(dim_param));
309 0 : data->fold_tight = isl_qpolynomial_fold_empty(data->type, dim_param);
310 0 : extract_coefficients(poly, dom, data);
311 :
312 0 : pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),
313 : data->fold);
314 0 : data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);
315 0 : pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, data->fold_tight);
316 0 : data->pwf_tight = isl_pw_qpolynomial_fold_fold(data->pwf_tight, pwf);
317 :
318 0 : isl_qpolynomial_free(poly);
319 0 : isl_cell_free(cell);
320 0 : for (i = 0; i < 1 + nvar; ++i)
321 0 : isl_qpolynomial_free(subs[i]);
322 0 : free(subs);
323 0 : return isl_stat_ok;
324 : error:
325 0 : isl_cell_free(cell);
326 0 : return isl_stat_error;
327 : }
328 :
329 : /* Base case of applying bernstein expansion.
330 : *
331 : * We compute the chamber decomposition of the parametric polytope "bset"
332 : * and then perform bernstein expansion on the parametric vertices
333 : * that are active on each chamber.
334 : */
335 0 : static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_base(
336 : __isl_take isl_basic_set *bset,
337 : __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
338 : {
339 : unsigned nvar;
340 : isl_space *dim;
341 : isl_pw_qpolynomial_fold *pwf;
342 : isl_vertices *vertices;
343 : int covers;
344 :
345 0 : nvar = isl_basic_set_dim(bset, isl_dim_set);
346 0 : if (nvar == 0) {
347 : isl_set *dom;
348 : isl_qpolynomial_fold *fold;
349 :
350 0 : fold = isl_qpolynomial_fold_alloc(data->type, poly);
351 0 : dom = isl_set_from_basic_set(bset);
352 0 : if (tight)
353 0 : *tight = 1;
354 0 : pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
355 0 : return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);
356 : }
357 :
358 0 : if (isl_qpolynomial_is_zero(poly)) {
359 : isl_set *dom;
360 : isl_qpolynomial_fold *fold;
361 0 : fold = isl_qpolynomial_fold_alloc(data->type, poly);
362 0 : dom = isl_set_from_basic_set(bset);
363 0 : pwf = isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
364 0 : if (tight)
365 0 : *tight = 1;
366 0 : return isl_pw_qpolynomial_fold_project_domain_on_params(pwf);
367 : }
368 :
369 0 : dim = isl_basic_set_get_space(bset);
370 0 : dim = isl_space_params(dim);
371 0 : dim = isl_space_from_domain(dim);
372 0 : dim = isl_space_add_dims(dim, isl_dim_set, 1);
373 0 : data->pwf = isl_pw_qpolynomial_fold_zero(isl_space_copy(dim), data->type);
374 0 : data->pwf_tight = isl_pw_qpolynomial_fold_zero(dim, data->type);
375 0 : data->poly = isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly));
376 0 : vertices = isl_basic_set_compute_vertices(bset);
377 0 : if (isl_vertices_foreach_disjoint_cell(vertices,
378 : &bernstein_coefficients_cell, data) < 0)
379 0 : data->pwf = isl_pw_qpolynomial_fold_free(data->pwf);
380 0 : isl_vertices_free(vertices);
381 0 : isl_qpolynomial_free(data->poly);
382 :
383 0 : isl_basic_set_free(bset);
384 0 : isl_qpolynomial_free(poly);
385 :
386 0 : covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);
387 0 : if (covers < 0)
388 0 : goto error;
389 :
390 0 : if (tight)
391 0 : *tight = covers;
392 :
393 0 : if (covers) {
394 0 : isl_pw_qpolynomial_fold_free(data->pwf);
395 0 : return data->pwf_tight;
396 : }
397 :
398 0 : data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);
399 :
400 0 : return data->pwf;
401 : error:
402 0 : isl_pw_qpolynomial_fold_free(data->pwf_tight);
403 0 : isl_pw_qpolynomial_fold_free(data->pwf);
404 0 : return NULL;
405 : }
406 :
407 : /* Apply bernstein expansion recursively by working in on len[i]
408 : * set variables at a time, with i ranging from n_group - 1 to 0.
409 : */
410 0 : static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_recursive(
411 : __isl_take isl_pw_qpolynomial *pwqp,
412 : int n_group, int *len, struct bernstein_data *data, int *tight)
413 : {
414 : int i;
415 : unsigned nparam;
416 : unsigned nvar;
417 : isl_pw_qpolynomial_fold *pwf;
418 :
419 0 : if (!pwqp)
420 0 : return NULL;
421 :
422 0 : nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
423 0 : nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_in);
424 :
425 0 : pwqp = isl_pw_qpolynomial_move_dims(pwqp, isl_dim_param, nparam,
426 0 : isl_dim_in, 0, nvar - len[n_group - 1]);
427 0 : pwf = isl_pw_qpolynomial_bound(pwqp, data->type, tight);
428 :
429 0 : for (i = n_group - 2; i >= 0; --i) {
430 0 : nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);
431 0 : pwf = isl_pw_qpolynomial_fold_move_dims(pwf, isl_dim_in, 0,
432 0 : isl_dim_param, nparam - len[i], len[i]);
433 0 : if (tight && !*tight)
434 0 : tight = NULL;
435 0 : pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);
436 : }
437 :
438 0 : return pwf;
439 : }
440 :
441 0 : static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_factors(
442 : __isl_take isl_basic_set *bset,
443 : __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
444 : {
445 : isl_factorizer *f;
446 : isl_set *set;
447 : isl_pw_qpolynomial *pwqp;
448 : isl_pw_qpolynomial_fold *pwf;
449 :
450 0 : f = isl_basic_set_factorizer(bset);
451 0 : if (!f)
452 0 : goto error;
453 0 : if (f->n_group == 0) {
454 0 : isl_factorizer_free(f);
455 0 : return bernstein_coefficients_base(bset, poly, data, tight);
456 : }
457 :
458 0 : set = isl_set_from_basic_set(bset);
459 0 : pwqp = isl_pw_qpolynomial_alloc(set, poly);
460 0 : pwqp = isl_pw_qpolynomial_morph_domain(pwqp, isl_morph_copy(f->morph));
461 :
462 0 : pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,
463 : tight);
464 :
465 0 : isl_factorizer_free(f);
466 :
467 0 : return pwf;
468 : error:
469 0 : isl_basic_set_free(bset);
470 0 : isl_qpolynomial_free(poly);
471 0 : return NULL;
472 : }
473 :
474 0 : static __isl_give isl_pw_qpolynomial_fold *bernstein_coefficients_full_recursive(
475 : __isl_take isl_basic_set *bset,
476 : __isl_take isl_qpolynomial *poly, struct bernstein_data *data, int *tight)
477 : {
478 : int i;
479 : int *len;
480 : unsigned nvar;
481 : isl_pw_qpolynomial_fold *pwf;
482 : isl_set *set;
483 : isl_pw_qpolynomial *pwqp;
484 :
485 0 : if (!bset || !poly)
486 : goto error;
487 :
488 0 : nvar = isl_basic_set_dim(bset, isl_dim_set);
489 :
490 0 : len = isl_alloc_array(bset->ctx, int, nvar);
491 0 : if (nvar && !len)
492 0 : goto error;
493 :
494 0 : for (i = 0; i < nvar; ++i)
495 0 : len[i] = 1;
496 :
497 0 : set = isl_set_from_basic_set(bset);
498 0 : pwqp = isl_pw_qpolynomial_alloc(set, poly);
499 :
500 0 : pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);
501 :
502 0 : free(len);
503 :
504 0 : return pwf;
505 : error:
506 0 : isl_basic_set_free(bset);
507 0 : isl_qpolynomial_free(poly);
508 0 : return NULL;
509 : }
510 :
511 : /* Compute a bound on the polynomial defined over the parametric polytope
512 : * using bernstein expansion and store the result
513 : * in bound->pwf and bound->pwf_tight.
514 : *
515 : * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if
516 : * the polytope can be factorized and apply bernstein expansion recursively
517 : * on the factors.
518 : * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply
519 : * bernstein expansion recursively on each dimension.
520 : * Otherwise, we apply bernstein expansion on the entire polytope.
521 : */
522 0 : isl_stat isl_qpolynomial_bound_on_domain_bernstein(
523 : __isl_take isl_basic_set *bset, __isl_take isl_qpolynomial *poly,
524 : struct isl_bound *bound)
525 : {
526 : struct bernstein_data data;
527 : isl_pw_qpolynomial_fold *pwf;
528 : unsigned nvar;
529 0 : int tight = 0;
530 0 : int *tp = bound->check_tight ? &tight : NULL;
531 :
532 0 : if (!bset || !poly)
533 : goto error;
534 :
535 0 : data.type = bound->type;
536 0 : data.check_tight = bound->check_tight;
537 :
538 0 : nvar = isl_basic_set_dim(bset, isl_dim_set);
539 :
540 0 : if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)
541 0 : pwf = bernstein_coefficients_factors(bset, poly, &data, tp);
542 0 : else if (nvar > 1 &&
543 0 : (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))
544 0 : pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);
545 : else
546 0 : pwf = bernstein_coefficients_base(bset, poly, &data, tp);
547 :
548 0 : if (tight)
549 0 : bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf);
550 : else
551 0 : bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf);
552 :
553 0 : return isl_stat_ok;
554 : error:
555 0 : isl_basic_set_free(bset);
556 0 : isl_qpolynomial_free(poly);
557 0 : return isl_stat_error;
558 : }
|
