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