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

source code of polly/lib/External/isl/isl_bernstein.c