nexmon – Blame information for rev 1

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