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_map_private.h>
  17 #include <isl/set.h>
  18 #include <isl/seq.h>
  19 #include <isl_morph.h>
  20 #include <isl_factorization.h>
  21 #include <isl_vertices_private.h>
  22 #include <isl_polynomial_private.h>
  23 #include <isl_bernstein.h>
  24
  25 struct bernstein_data {
  26         enum isl_fold type;
  27         isl_qpolynomial *poly;
  28         int check_tight;
  29
  30         isl_cell *cell;
  31
  32         isl_qpolynomial_fold *fold;
  33         isl_qpolynomial_fold *fold_tight;
  34         isl_pw_qpolynomial_fold *pwf;
  35         isl_pw_qpolynomial_fold *pwf_tight;
  36 };
  37
  38 static int vertex_is_integral(__isl_keep isl_basic_set *vertex)
  39 {
  40         unsigned nvar;
  41         unsigned nparam;
  42         int i;
  43
  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) {
  47                 int r = nvar - 1 - i;
  48                 if (!isl_int_is_one(vertex->eq[r][1 + nparam + i]) &&
  49                     !isl_int_is_negone(vertex->eq[r][1 + nparam + i]))
  50                         return 0;
  51         }
  52
  53         return 1;
  54 }
  55
  56 static __isl_give isl_qpolynomial *vertex_coordinate(
  57         __isl_keep isl_basic_set *vertex, int i, __isl_take isl_dim *dim)
  58 {
  59         unsigned nvar;
  60         unsigned nparam;
  61         int r;
  62         isl_int denom;
  63         isl_qpolynomial *v;
  64
  65         nvar = isl_basic_set_dim(vertex, isl_dim_set);
  66         nparam = isl_basic_set_dim(vertex, isl_dim_param);
  67         r = nvar - 1 - i;
  68
  69         isl_int_init(denom);
  70         isl_int_set(denom, vertex->eq[r][1 + nparam + i]);
  71         isl_assert(vertex->ctx, !isl_int_is_zero(denom), goto error);
  72
  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));
  76         else
  77                 isl_int_neg(denom, denom);
  78
  79         v = isl_qpolynomial_from_affine(dim, vertex->eq[r], denom);
  80         isl_int_clear(denom);
  81
  82         return v;
  83 error:
  84         isl_dim_free(dim);
  85         isl_int_clear(denom);
  86         return NULL;
  87 }
  88
  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.
  92  */
  93 static int is_tight(int *k, int n, int d, isl_cell *cell)
  94 {
  95         int i, j;
  96
  97         for (i = 0; i < n; ++i) {
  98                 int v;
  99                 if (k[i] != d) {
 100                         if (k[i])
 101                                 return 0;
 102                         continue;
 103                 }
 104                 v = cell->ids[n - 1 - i];
 105                 return vertex_is_integral(cell->vertices->v[v].vertex);
 106         }
 107
 108         return 0;
 109 }
 110
 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)
 113 {
 114         isl_qpolynomial_fold *fold;
 115
 116         fold = isl_qpolynomial_fold_alloc(data->type, b);
 117
 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);
 121         else
 122                 data->fold = isl_qpolynomial_fold_fold_on_domain(dom,
 123                                                         data->fold, fold);
 124 }
 125
 126 /* Extract the coefficients of the Bernstein base polynomials and store
 127  * them in data->fold and data->fold_tight.
 128  *
 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]!
 132  *
 133  * c[i] contains the coefficient of the selected powers of the first i+1 vars.
 134  * multinom[i] contains the partial multinomial coefficient.
 135  */
 136 static void extract_coefficients(isl_qpolynomial *poly,
 137         __isl_keep isl_set *dom, struct bernstein_data *data)
 138 {
 139         int i;
 140         int d;
 141         int n;
 142         isl_ctx *ctx;
 143         isl_qpolynomial **c = NULL;
 144         int *k = NULL;
 145         int *left = NULL;
 146         isl_vec *multinom = NULL;
 147
 148         if (!poly)
 149                 return;
 150
 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);
 155
 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)
 161                 goto error;
 162
 163         isl_int_set_si(multinom->el[0], 1);
 164         for (k[0] = d; k[0] >= 0; --k[0]) {
 165                 int i = 1;
 166                 isl_qpolynomial_free(c[0]);
 167                 c[0] = isl_qpolynomial_coeff(poly, isl_dim_set, n - 1, k[0]);
 168                 left[0] = d - k[0];
 169                 k[1] = -1;
 170                 isl_int_set(multinom->el[1], multinom->el[0]);
 171                 while (i > 0) {
 172                         if (i == n - 1) {
 173                                 int j;
 174                                 isl_dim *dim;
 175                                 isl_qpolynomial *b;
 176                                 isl_qpolynomial *f;
 177                                 for (j = 2; j <= left[i - 1]; ++j)
 178                                         isl_int_divexact_ui(multinom->el[i],
 179                                                 multinom->el[i], j);
 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,
 183                                                                 0, n);
 184                                 dim = isl_qpolynomial_get_dim(b);
 185                                 f = isl_qpolynomial_rat_cst(dim, ctx->one,
 186                                         multinom->el[i]);
 187                                 b = isl_qpolynomial_mul(b, f);
 188                                 k[n - 1] = left[n - 2];
 189                                 add_fold(b, dom, k, n, d, data);
 190                                 --i;
 191                                 continue;
 192                         }
 193                         if (k[i] >= left[i - 1]) {
 194                                 --i;
 195                                 continue;
 196                         }
 197                         ++k[i];
 198                         if (k[i])
 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,
 203                                         n - 1 - i, k[i]);
 204                         left[i] = left[i - 1] - k[i];
 205                         k[i + 1] = -1;
 206                         isl_int_set(multinom->el[i + 1], multinom->el[i]);
 207                         ++i;
 208                 }
 209                 isl_int_mul_ui(multinom->el[0], multinom->el[0], k[0]);
 210         }
 211
 212         for (i = 0; i < n; ++i)
 213                 isl_qpolynomial_free(c[i]);
 214
 215         isl_vec_free(multinom);
 216         free(left);
 217         free(k);
 218         free(c);
 219         return;
 220 error:
 221         isl_vec_free(multinom);
 222         free(left);
 223         free(k);
 224         if (c)
 225                 for (i = 0; i < n; ++i)
 226                         isl_qpolynomial_free(c[i]);
 227         free(c);
 228         return;
 229 }
 230
 231 /* Perform bernstein expansion on the parametric vertices that are active
 232  * on "cell".
 233  *
 234  * data->poly has been homogenized in the calling function.
 235  *
 236  * We plug in the barycentric coordinates for the set variables
 237  *
 238  *              \vec x = \sum_i \alpha_i v_i(\vec p)
 239  *
 240  * and the constant "1 = \sum_i \alpha_i" for the homogeneous dimension.
 241  * Next, we extract the coefficients of the Bernstein base polynomials.
 242  */
 243 static int bernstein_coefficients_cell(__isl_take isl_cell *cell, void *user)
 244 {
 245         int i, j;
 246         struct bernstein_data *data = (struct bernstein_data *)user;
 247         isl_dim *dim_param;
 248         isl_dim *dim_dst;
 249         isl_qpolynomial *poly = data->poly;
 250         unsigned nvar;
 251         int n_vertices;
 252         isl_qpolynomial **subs;
 253         isl_pw_qpolynomial_fold *pwf;
 254         isl_set *dom;
 255         isl_ctx *ctx;
 256
 257         nvar = isl_qpolynomial_dim(poly, isl_dim_set) - 1;
 258         n_vertices = cell->n_vertices;
 259
 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);
 264
 265         subs = isl_alloc_array(data->poly->dim->ctx, isl_qpolynomial *,
 266                                 1 + nvar);
 267         if (!subs)
 268                 goto error;
 269
 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);
 273
 274         for (i = 0; i < 1 + nvar; ++i)
 275                 subs[i] = isl_qpolynomial_zero(isl_dim_copy(dim_dst));
 276
 277         for (i = 0; i < n_vertices; ++i) {
 278                 isl_qpolynomial *c;
 279                 c = isl_qpolynomial_var(isl_dim_copy(dim_dst), isl_dim_set,
 280                                         1 + nvar + i);
 281                 for (j = 0; j < nvar; ++j) {
 282                         int k = cell->ids[i];
 283                         isl_qpolynomial *v;
 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);
 290                 }
 291                 subs[0] = isl_qpolynomial_add(subs[0], c);
 292         }
 293         isl_dim_free(dim_dst);
 294
 295         poly = isl_qpolynomial_copy(poly);
 296
 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);
 300
 301         data->cell = cell;
 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);
 306
 307         pwf = isl_pw_qpolynomial_fold_alloc(data->type, isl_set_copy(dom),
 308                                             data->fold);
 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);
 312
 313         isl_qpolynomial_free(poly);
 314         isl_cell_free(cell);
 315         for (i = 0; i < 1 + nvar; ++i)
 316                 isl_qpolynomial_free(subs[i]);
 317         free(subs);
 318         return 0;
 319 error:
 320         isl_cell_free(cell);
 321         return -1;
 322 }
 323
 324 /* Base case of applying bernstein expansion.
 325  *
 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.
 329  */
 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)
 333 {
 334         unsigned nvar;
 335         isl_dim *dim;
 336         isl_pw_qpolynomial_fold *pwf;
 337         isl_vertices *vertices;
 338         int covers;
 339
 340         nvar = isl_basic_set_dim(bset, isl_dim_set);
 341         if (nvar == 0) {
 342                 isl_set *dom;
 343                 isl_qpolynomial_fold *fold;
 344                 fold = isl_qpolynomial_fold_alloc(data->type, poly);
 345                 dom = isl_set_from_basic_set(bset);
 346                 if (tight)
 347                         *tight = 1;
 348                 return isl_pw_qpolynomial_fold_alloc(data->type, dom, fold);
 349         }
 350
 351         if (isl_qpolynomial_is_zero(poly)) {
 352                 isl_set *dom;
 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);
 357                 if (tight)
 358                         *tight = 1;
 359                 return isl_pw_qpolynomial_fold_drop_dims(pwf,
 360                                                             isl_dim_set, 0, nvar);
 361         }
 362
 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);
 373
 374         isl_basic_set_free(bset);
 375         isl_qpolynomial_free(poly);
 376
 377         covers = isl_pw_qpolynomial_fold_covers(data->pwf_tight, data->pwf);
 378         if (covers < 0)
 379                 goto error;
 380
 381         if (tight)
 382                 *tight = covers;
 383
 384         if (covers) {
 385                 isl_pw_qpolynomial_fold_free(data->pwf);
 386                 return data->pwf_tight;
 387         }
 388
 389         data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, data->pwf_tight);
 390
 391         return data->pwf;
 392 error:
 393         isl_pw_qpolynomial_fold_free(data->pwf_tight);
 394         isl_pw_qpolynomial_fold_free(data->pwf);
 395         return NULL;
 396 }
 397
 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.
 400  */
 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)
 404 {
 405         int i;
 406         unsigned nparam;
 407         unsigned nvar;
 408         isl_pw_qpolynomial_fold *pwf;
 409
 410         if (!pwqp)
 411                 return NULL;
 412
 413         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 414         nvar = isl_pw_qpolynomial_dim(pwqp, isl_dim_set);
 415
 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);
 419
 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)
 425                         tight = NULL;
 426                 pwf = isl_pw_qpolynomial_fold_bound(pwf, tight);
 427         }
 428
 429         return pwf;
 430 }
 431
 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)
 435 {
 436         isl_factorizer *f;
 437         isl_set *set;
 438         isl_pw_qpolynomial *pwqp;
 439         isl_pw_qpolynomial_fold *pwf;
 440
 441         f = isl_basic_set_factorizer(bset);
 442         if (!f)
 443                 goto error;
 444         if (f->n_group == 0) {
 445                 isl_factorizer_free(f);
 446                 return  bernstein_coefficients_base(bset, poly, data, tight);
 447         }
 448
 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));
 452
 453         pwf = bernstein_coefficients_recursive(pwqp, f->n_group, f->len, data,
 454                                                 tight);
 455
 456         isl_factorizer_free(f);
 457
 458         return pwf;
 459 error:
 460         isl_basic_set_free(bset);
 461         isl_qpolynomial_free(poly);
 462         return NULL;
 463 }
 464
 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)
 468 {
 469         int i;
 470         int *len;
 471         unsigned nvar;
 472         isl_pw_qpolynomial_fold *pwf;
 473         isl_set *set;
 474         isl_pw_qpolynomial *pwqp;
 475
 476         if (!bset || !poly)
 477                 goto error;
 478
 479         nvar = isl_basic_set_dim(bset, isl_dim_set);
 480
 481         len = isl_alloc_array(bset->ctx, int, nvar);
 482         if (!len)
 483                 goto error;
 484
 485         for (i = 0; i < nvar; ++i)
 486                 len[i] = 1;
 487
 488         set = isl_set_from_basic_set(bset);
 489         pwqp = isl_pw_qpolynomial_alloc(set, poly);
 490
 491         pwf = bernstein_coefficients_recursive(pwqp, nvar, len, data, tight);
 492
 493         free(len);
 494
 495         return pwf;
 496 error:
 497         isl_basic_set_free(bset);
 498         isl_qpolynomial_free(poly);
 499         return NULL;
 500 }
 501
 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.
 505  *
 506  * If bernstein_recurse is set to ISL_BERNSTEIN_FACTORS, we check if
 507  * the polytope can be factorized and apply bernstein expansion recursively
 508  * on the factors.
 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.
 512  */
 513 int isl_qpolynomial_bound_on_domain_bernstein(__isl_take isl_basic_set *bset,
 514         __isl_take isl_qpolynomial *poly, struct isl_bound *bound)
 515 {
 516         struct bernstein_data data;
 517         isl_pw_qpolynomial_fold *pwf;
 518         unsigned nvar;
 519         int tight = 0;
 520         int *tp = bound->check_tight ? &tight : NULL;
 521
 522         if (!bset || !poly)
 523                 goto error;
 524
 525         data.type = bound->type;
 526         data.check_tight = bound->check_tight;
 527
 528         nvar = isl_basic_set_dim(bset, isl_dim_set);
 529
 530         if (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_FACTORS)
 531                 pwf = bernstein_coefficients_factors(bset, poly, &data, tp);
 532         else if (nvar > 1 &&
 533             (bset->ctx->opt->bernstein_recurse & ISL_BERNSTEIN_INTERVALS))
 534                 pwf = bernstein_coefficients_full_recursive(bset, poly, &data, tp);
 535         else
 536                 pwf = bernstein_coefficients_base(bset, poly, &data, tp);
 537
 538         if (tight)
 539                 bound->pwf_tight = isl_pw_qpolynomial_fold_fold(bound->pwf_tight, pwf);
 540         else
 541                 bound->pwf = isl_pw_qpolynomial_fold_fold(bound->pwf, pwf);
 542
 543         return 0;
 544 error:
 545         isl_basic_set_free(bset);
 546         isl_qpolynomial_free(poly);
 547         return -1;
 548 }