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
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>
23 #include <isl_map_private.h>
25 struct bernstein_data
{
27 isl_qpolynomial
*poly
;
32 isl_qpolynomial_fold
*fold
;
33 isl_qpolynomial_fold
*fold_tight
;
34 isl_pw_qpolynomial_fold
*pwf
;
35 isl_pw_qpolynomial_fold
*pwf_tight
;
38 static int vertex_is_integral(__isl_keep isl_basic_set
*vertex
)
44 nvar
= isl_basic_set_dim(vertex
, isl_dim_set
);
45 nparam
= isl_basic_set_dim(vertex
, isl_dim_param
);
46 for (i
= 0; i
< nvar
; ++i
) {
48 if (!isl_int_is_one(vertex
->eq
[r
][1 + nparam
+ i
]) &&
49 !isl_int_is_negone(vertex
->eq
[r
][1 + nparam
+ i
]))
56 static __isl_give isl_qpolynomial
*vertex_coordinate(
57 __isl_keep isl_basic_set
*vertex
, int i
, __isl_take isl_dim
*dim
)
65 nvar
= isl_basic_set_dim(vertex
, isl_dim_set
);
66 nparam
= isl_basic_set_dim(vertex
, isl_dim_param
);
70 isl_int_set(denom
, vertex
->eq
[r
][1 + nparam
+ i
]);
71 isl_assert(vertex
->ctx
, !isl_int_is_zero(denom
), goto error
);
73 if (isl_int_is_pos(denom
))
74 isl_seq_neg(vertex
->eq
[r
], vertex
->eq
[r
],
75 1 + isl_basic_set_total_dim(vertex
));
77 isl_int_neg(denom
, denom
);
79 v
= isl_qpolynomial_from_affine(dim
, vertex
->eq
[r
], denom
);
89 /* Check whether the bound associated to the selection "k" is tight,
90 * which is the case if we select exactly one vertex and if that vertex
91 * is integral for all values of the parameters.
93 static int is_tight(int *k
, int n
, int d
, isl_cell
*cell
)
97 for (i
= 0; i
< n
; ++i
) {
104 v
= cell
->ids
[n
- 1 - i
];
105 return vertex_is_integral(cell
->vertices
->v
[v
].vertex
);
111 static void add_fold(__isl_take isl_qpolynomial
*b
, __isl_keep isl_set
*dom
,
112 int *k
, int n
, int d
, struct bernstein_data
*data
)
114 isl_qpolynomial_fold
*fold
;
116 fold
= isl_qpolynomial_fold_alloc(data
->type
, b
);
118 if (data
->check_tight
&& is_tight(k
, n
, d
, data
->cell
))
119 data
->fold_tight
= isl_qpolynomial_fold_fold_on_domain(dom
,
120 data
->fold_tight
, fold
);
122 data
->fold
= isl_qpolynomial_fold_fold_on_domain(dom
,
126 /* Extract the coefficients of the Bernstein base polynomials and store
127 * them in data->fold and data->fold_tight.
129 * In particular, the coefficient of each monomial
130 * of multi-degree (k[0], k[1], ..., k[n-1]) is divided by the corresponding
131 * multinomial coefficient d!/k[0]! k[1]! ... k[n-1]!
133 * c[i] contains the coefficient of the selected powers of the first i+1 vars.
134 * multinom[i] contains the partial multinomial coefficient.
136 static void extract_coefficients(isl_qpolynomial
*poly
,
137 __isl_keep isl_set
*dom
, struct bernstein_data
*data
)
143 isl_qpolynomial
**c
= NULL
;
146 isl_vec
*multinom
= NULL
;
151 ctx
= isl_qpolynomial_get_ctx(poly
);
152 n
= isl_qpolynomial_dim(poly
, isl_dim_set
);
153 d
= isl_qpolynomial_degree(poly
);
154 isl_assert(ctx
, n
>= 2, return);
156 c
= isl_calloc_array(ctx
, isl_qpolynomial
*, n
);
157 k
= isl_alloc_array(ctx
, int, n
);
158 left
= isl_alloc_array(ctx
, int, n
);
159 multinom
= isl_vec_alloc(ctx
, n
);
160 if (!c
|| !k
|| !left
|| !multinom
)
163 isl_int_set_si(multinom
->el
[0], 1);
164 for (k
[0] = d
; k
[0] >= 0; --k
[0]) {
166 isl_qpolynomial_free(c
[0]);
167 c
[0] = isl_qpolynomial_coeff(poly
, isl_dim_set
, n
- 1, k
[0]);
170 isl_int_set(multinom
->el
[1], multinom
->el
[0]);
177 for (j
= 2; j
<= left
[i
- 1]; ++j
)
178 isl_int_divexact_ui(multinom
->el
[i
],
180 b
= isl_qpolynomial_coeff(c
[i
- 1], isl_dim_set
,
181 n
- 1 - i
, left
[i
- 1]);
182 b
= isl_qpolynomial_drop_dims(b
, isl_dim_set
,
184 dim
= isl_qpolynomial_get_dim(b
);
185 f
= isl_qpolynomial_rat_cst(dim
, ctx
->one
,
187 b
= isl_qpolynomial_mul(b
, f
);
188 k
[n
- 1] = left
[n
- 2];
189 add_fold(b
, dom
, k
, n
, d
, data
);
193 if (k
[i
] >= left
[i
- 1]) {
199 isl_int_divexact_ui(multinom
->el
[i
],
200 multinom
->el
[i
], k
[i
]);
201 isl_qpolynomial_free(c
[i
]);
202 c
[i
] = isl_qpolynomial_coeff(c
[i
- 1], isl_dim_set
,
204 left
[i
] = left
[i
- 1] - k
[i
];
206 isl_int_set(multinom
->el
[i
+ 1], multinom
->el
[i
]);
209 isl_int_mul_ui(multinom
->el
[0], multinom
->el
[0], k
[0]);
212 for (i
= 0; i
< n
; ++i
)
213 isl_qpolynomial_free(c
[i
]);
215 isl_vec_free(multinom
);
221 isl_vec_free(multinom
);
225 for (i
= 0; i
< n
; ++i
)
226 isl_qpolynomial_free(c
[i
]);
231 /* Perform bernstein expansion on the parametric vertices that are active
234 * data->poly has been homogenized in the calling function.
236 * We plug in the barycentric coordinates for the set variables
238 * \vec x = \sum_i \alpha_i v_i(\vec p)
240 * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.
241 * Next, we extract the coefficients of the Bernstein base polynomials.
243 static int bernstein_coefficients_cell(__isl_take isl_cell
*cell
, void *user
)
246 struct bernstein_data
*data
= (struct bernstein_data
*)user
;
249 isl_qpolynomial
*poly
= data
->poly
;
252 isl_qpolynomial
**subs
;
253 isl_pw_qpolynomial_fold
*pwf
;
257 nvar
= isl_qpolynomial_dim(poly
, isl_dim_set
) - 1;
258 n_vertices
= cell
->n_vertices
;
260 ctx
= isl_qpolynomial_get_ctx(poly
);
261 if (n_vertices
> nvar
+ 1 && ctx
->opt
->bernstein_triangulate
)
262 return isl_cell_foreach_simplex(cell
,
263 &bernstein_coefficients_cell
, user
);
265 subs
= isl_alloc_array(data
->poly
->dim
->ctx
, isl_qpolynomial
*,
270 dim_param
= isl_basic_set_get_dim(cell
->dom
);
271 dim_dst
= isl_qpolynomial_get_dim(poly
);
272 dim_dst
= isl_dim_add(dim_dst
, isl_dim_set
, n_vertices
);
274 for (i
= 0; i
< 1 + nvar
; ++i
)
275 subs
[i
] = isl_qpolynomial_zero(isl_dim_copy(dim_dst
));
277 for (i
= 0; i
< n_vertices
; ++i
) {
279 c
= isl_qpolynomial_var(isl_dim_copy(dim_dst
), isl_dim_set
,
281 for (j
= 0; j
< nvar
; ++j
) {
282 int k
= cell
->ids
[i
];
284 v
= vertex_coordinate(cell
->vertices
->v
[k
].vertex
, j
,
285 isl_dim_copy(dim_param
));
286 v
= isl_qpolynomial_add_dims(v
, isl_dim_set
,
287 1 + nvar
+ n_vertices
);
288 v
= isl_qpolynomial_mul(v
, isl_qpolynomial_copy(c
));
289 subs
[1 + j
] = isl_qpolynomial_add(subs
[1 + j
], v
);
291 subs
[0] = isl_qpolynomial_add(subs
[0], c
);
293 isl_dim_free(dim_dst
);
295 poly
= isl_qpolynomial_copy(poly
);
297 poly
= isl_qpolynomial_add_dims(poly
, isl_dim_set
, n_vertices
);
298 poly
= isl_qpolynomial_substitute(poly
, isl_dim_set
, 0, 1 + nvar
, subs
);
299 poly
= isl_qpolynomial_drop_dims(poly
, isl_dim_set
, 0, 1 + nvar
);
302 dom
= isl_set_from_basic_set(isl_basic_set_copy(cell
->dom
));
303 data
->fold
= isl_qpolynomial_fold_empty(data
->type
, isl_dim_copy(dim_param
));
304 data
->fold_tight
= isl_qpolynomial_fold_empty(data
->type
, dim_param
);
305 extract_coefficients(poly
, dom
, data
);
307 pwf
= isl_pw_qpolynomial_fold_alloc(data
->type
, isl_set_copy(dom
),
309 data
->pwf
= isl_pw_qpolynomial_fold_fold(data
->pwf
, pwf
);
310 pwf
= isl_pw_qpolynomial_fold_alloc(data
->type
, dom
, data
->fold_tight
);
311 data
->pwf_tight
= isl_pw_qpolynomial_fold_fold(data
->pwf_tight
, pwf
);
313 isl_qpolynomial_free(poly
);
315 for (i
= 0; i
< 1 + nvar
; ++i
)
316 isl_qpolynomial_free(subs
[i
]);
324 /* Base case of applying bernstein expansion.
326 * We compute the chamber decomposition of the parametric polytope "bset"
327 * and then perform bernstein expansion on the parametric vertices
328 * that are active on each chamber.
330 static __isl_give isl_pw_qpolynomial_fold
*bernstein_coefficients_base(
331 __isl_take isl_basic_set
*bset
,
332 __isl_take isl_qpolynomial
*poly
, struct bernstein_data
*data
, int *tight
)
336 isl_pw_qpolynomial_fold
*pwf
;
337 isl_vertices
*vertices
;
340 nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
343 isl_qpolynomial_fold
*fold
;
344 fold
= isl_qpolynomial_fold_alloc(data
->type
, poly
);
345 dom
= isl_set_from_basic_set(bset
);
348 return isl_pw_qpolynomial_fold_alloc(data
->type
, dom
, fold
);
351 if (isl_qpolynomial_is_zero(poly
)) {
353 isl_qpolynomial_fold
*fold
;
354 fold
= isl_qpolynomial_fold_alloc(data
->type
, poly
);
355 dom
= isl_set_from_basic_set(bset
);
356 pwf
= isl_pw_qpolynomial_fold_alloc(data
->type
, dom
, fold
);
359 return isl_pw_qpolynomial_fold_drop_dims(pwf
,
360 isl_dim_set
, 0, nvar
);
363 dim
= isl_basic_set_get_dim(bset
);
364 dim
= isl_dim_drop(dim
, isl_dim_set
, 0, nvar
);
365 data
->pwf
= isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim
), data
->type
);
366 data
->pwf_tight
= isl_pw_qpolynomial_fold_zero(dim
, data
->type
);
367 data
->poly
= isl_qpolynomial_homogenize(isl_qpolynomial_copy(poly
));
368 vertices
= isl_basic_set_compute_vertices(bset
);
369 isl_vertices_foreach_disjoint_cell(vertices
,
370 &bernstein_coefficients_cell
, data
);
371 isl_vertices_free(vertices
);
372 isl_qpolynomial_free(data
->poly
);
374 isl_basic_set_free(bset
);
375 isl_qpolynomial_free(poly
);
377 covers
= isl_pw_qpolynomial_fold_covers(data
->pwf_tight
, data
->pwf
);
385 isl_pw_qpolynomial_fold_free(data
->pwf
);
386 return data
->pwf_tight
;
389 data
->pwf
= isl_pw_qpolynomial_fold_fold(data
->pwf
, data
->pwf_tight
);
393 isl_pw_qpolynomial_fold_free(data
->pwf_tight
);
394 isl_pw_qpolynomial_fold_free(data
->pwf
);
398 /* Apply bernstein expansion recursively by working in on len[i]
399 * set variables at a time, with i ranging from n_group - 1 to 0.
401 static __isl_give isl_pw_qpolynomial_fold
*bernstein_coefficients_recursive(
402 __isl_take isl_pw_qpolynomial
*pwqp
,
403 int n_group
, int *len
, struct bernstein_data
*data
, int *tight
)
408 isl_pw_qpolynomial_fold
*pwf
;
413 nparam
= isl_pw_qpolynomial_dim(pwqp
, isl_dim_param
);
414 nvar
= isl_pw_qpolynomial_dim(pwqp
, isl_dim_set
);
416 pwqp
= isl_pw_qpolynomial_move_dims(pwqp
, isl_dim_param
, nparam
,
417 isl_dim_set
, 0, nvar
- len
[n_group
- 1]);
418 pwf
= isl_pw_qpolynomial_bound(pwqp
, data
->type
, tight
);
420 for (i
= n_group
- 2; i
>= 0; --i
) {
421 nparam
= isl_pw_qpolynomial_fold_dim(pwf
, isl_dim_param
);
422 pwf
= isl_pw_qpolynomial_fold_move_dims(pwf
, isl_dim_set
, 0,
423 isl_dim_param
, nparam
- len
[i
], len
[i
]);
424 if (tight
&& !*tight
)
426 pwf
= isl_pw_qpolynomial_fold_bound(pwf
, tight
);
432 static __isl_give isl_pw_qpolynomial_fold
*bernstein_coefficients_factors(
433 __isl_take isl_basic_set
*bset
,
434 __isl_take isl_qpolynomial
*poly
, struct bernstein_data
*data
, int *tight
)
438 isl_pw_qpolynomial
*pwqp
;
439 isl_pw_qpolynomial_fold
*pwf
;
441 f
= isl_basic_set_factorizer(bset
);
444 if (f
->n_group
== 0) {
445 isl_factorizer_free(f
);
446 return bernstein_coefficients_base(bset
, poly
, data
, tight
);
449 set
= isl_set_from_basic_set(bset
);
450 pwqp
= isl_pw_qpolynomial_alloc(set
, poly
);
451 pwqp
= isl_pw_qpolynomial_morph(pwqp
, isl_morph_copy(f
->morph
));
453 pwf
= bernstein_coefficients_recursive(pwqp
, f
->n_group
, f
->len
, data
,
456 isl_factorizer_free(f
);
460 isl_basic_set_free(bset
);
461 isl_qpolynomial_free(poly
);
465 static __isl_give isl_pw_qpolynomial_fold
*bernstein_coefficients_full_recursive(
466 __isl_take isl_basic_set
*bset
,
467 __isl_take isl_qpolynomial
*poly
, struct bernstein_data
*data
, int *tight
)
472 isl_pw_qpolynomial_fold
*pwf
;
474 isl_pw_qpolynomial
*pwqp
;
479 nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
481 len
= isl_alloc_array(bset
->ctx
, int, nvar
);
485 for (i
= 0; i
< nvar
; ++i
)
488 set
= isl_set_from_basic_set(bset
);
489 pwqp
= isl_pw_qpolynomial_alloc(set
, poly
);
491 pwf
= bernstein_coefficients_recursive(pwqp
, nvar
, len
, data
, tight
);
497 isl_basic_set_free(bset
);
498 isl_qpolynomial_free(poly
);
502 /* Compute a bound on the polynomial defined over the parametric polytope
503 * using bernstein expansion and store the result
504 * in bound->pwf and bound->pwf_tight.
506 * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if
507 * the polytope can be factorized and apply bernstein expansion recursively
509 * If bernstein_recurse is set to ISL_BERNSTEIN_INTERVALS, we apply
510 * bernstein expansion recursively on each dimension.
511 * Otherwise, we apply bernstein expansion on the entire polytope.
513 int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set
*bset
,
514 __isl_take isl_qpolynomial
*poly
, struct isl_bound
*bound
)
516 struct bernstein_data data
;
517 isl_pw_qpolynomial_fold
*pwf
;
520 int *tp
= bound
->check_tight
? &tight
: NULL
;
525 data
.type
= bound
->type
;
526 data
.check_tight
= bound
->check_tight
;
528 nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
530 if (bset
->ctx
->opt
->bernstein_recurse
& ISL_BERNSTEIN_FACTORS
)
531 pwf
= bernstein_coefficients_factors(bset
, poly
, &data
, tp
);
533 (bset
->ctx
->opt
->bernstein_recurse
& ISL_BERNSTEIN_INTERVALS
))
534 pwf
= bernstein_coefficients_full_recursive(bset
, poly
, &data
, tp
);
536 pwf
= bernstein_coefficients_base(bset
, poly
, &data
, tp
);
539 bound
->pwf_tight
= isl_pw_qpolynomial_fold_fold(bound
->pwf_tight
, pwf
);
541 bound
->pwf
= isl_pw_qpolynomial_fold_fold(bound
->pwf
, pwf
);
545 isl_basic_set_free(bset
);
546 isl_qpolynomial_free(poly
);