isl_bernstein.c

   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 }