bernstein.cc

   1 #include <assert.h>
   2 #include <vector>
   3 #include <bernstein/bernstein.h>
   4 #include <bernstein/piecewise_lst.h>
   5 #include <barvinok/barvinok.h>
   6 #include <barvinok/util.h>
   7 #include <barvinok/bernstein.h>
   8 #include <barvinok/options.h>
   9
  10 using namespace GiNaC;
  11 using namespace bernstein;
  12
  13 using std::pair;
  14 using std::vector;
  15 using std::cerr;
  16 using std::endl;
  17
  18 namespace barvinok {
  19
  20 ex evalue2ex(evalue *e, const exvector& vars)
  21 {
  22     if (value_notzero_p(e->d))
  23         return value2numeric(e->x.n)/value2numeric(e->d);
  24     if (e->x.p->type != polynomial)
  25         return fail();
  26     ex poly = 0;
  27     for (int i = e->x.p->size-1; i >= 0; --i) {
  28         poly *= vars[e->x.p->pos-1];
  29         ex t = evalue2ex(&e->x.p->arr[i], vars);
  30         if (is_exactly_a<fail>(t))
  31             return t;
  32         poly += t;
  33     }
  34     return poly;
  35 }
  36
  37 static int type_offset(enode *p)
  38 {
  39    return p->type == fractional ? 1 :
  40           p->type == flooring ? 1 : 0;
  41 }
  42
  43 typedef pair<bool, const evalue *> typed_evalue;
  44
  45 static ex evalue2ex_add_var(evalue *e, exvector& extravar,
  46                             vector<typed_evalue>& expr, bool is_fract)
  47 {
  48     ex base_var = 0;
  49
  50     for (int i = 0; i < expr.size(); ++i) {
  51         if (is_fract == expr[i].first && eequal(e, expr[i].second)) {
  52             base_var = extravar[i];
  53             break;
  54         }
  55     }
  56     if (base_var != 0)
  57         return base_var;
  58
  59     char name[20];
  60     snprintf(name, sizeof(name), "f%c%d", is_fract ? 'r' : 'l', expr.size());
  61     extravar.push_back(base_var = symbol(name));
  62     expr.push_back(typed_evalue(is_fract, e));
  63
  64     return base_var;
  65 }
  66
  67 /* For the argument e=(f/d) of a fractional, return (d-1)/d times
  68  * a variable in [0,1] (see setup_constraints).
  69  */
  70 static ex evalue2ex_get_fract(evalue *e, exvector& extravar,
  71                                 vector<typed_evalue>& expr)
  72 {
  73     ex f;
  74     Value d;
  75     ex den;
  76     value_init(d);
  77     value_set_si(d, 1);
  78     evalue_denom(e, &d);
  79     den = value2numeric(d);
  80     value_clear(d);
  81     f = (den-1)/den;
  82
  83     ex base_var = evalue2ex_add_var(e, extravar, expr, true);
  84     base_var *= f;
  85     return base_var;
  86 }
  87
  88 static ex evalue2ex_r(const evalue *e, const exvector& vars,
  89                       exvector& extravar, vector<typed_evalue>& expr,
  90                       Vector *coset)
  91 {
  92     if (value_notzero_p(e->d))
  93         return value2numeric(e->x.n)/value2numeric(e->d);
  94     ex base_var = 0;
  95     ex poly = 0;
  96
  97     switch (e->x.p->type) {
  98     case polynomial:
  99         base_var = vars[e->x.p->pos-1];
 100         break;
 101     case flooring:
 102         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, false);
 103         break;
 104     case fractional:
 105         base_var = evalue2ex_get_fract(&e->x.p->arr[0], extravar, expr);
 106         break;
 107     case periodic:
 108         assert(coset);
 109         return evalue2ex_r(&e->x.p->arr[VALUE_TO_INT(coset->p[e->x.p->pos-1])],
 110                            vars, extravar, expr, coset);
 111     default:
 112         return fail();
 113     }
 114
 115     int offset = type_offset(e->x.p);
 116     for (int i = e->x.p->size-1; i >= offset; --i) {
 117         poly *= base_var;
 118         ex t = evalue2ex_r(&e->x.p->arr[i], vars, extravar, expr, coset);
 119         if (is_exactly_a<fail>(t))
 120             return t;
 121         poly += t;
 122     }
 123     return poly;
 124 }
 125
 126 /* For each t = floor(e/d), set up two constraints
 127  *
 128  *       e - d t       >= 0
 129  *      -e + d t + d-1 >= 0
 130  *
 131  * e is assumed to be an affine expression.
 132  *
 133  * For each t = fract(e/d), set up two constraints
 134  *
 135  *          -d t + d-1 >= 0
 136  *             t       >= 0
 137  */
 138 static Matrix *setup_constraints(const vector<typed_evalue> expr, int nvar)
 139 {
 140     int extra = expr.size();
 141     if (!extra)
 142         return NULL;
 143     Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
 144     for (int i = 0; i < extra; ++i) {
 145         if (expr[i].first) {
 146             value_set_si(M->p[2*i][0], 1);
 147             value_set_si(M->p[2*i][1+i], -1);
 148             value_set_si(M->p[2*i][1+extra+nvar], 1);
 149             value_set_si(M->p[2*i+1][0], 1);
 150             value_set_si(M->p[2*i+1][1+i], 1);
 151         } else {
 152             Value *d = &M->p[2*i][1+i];
 153             evalue_extract_affine(expr[i].second, M->p[2*i]+1+extra,
 154                                   M->p[2*i]+1+extra+nvar, d);
 155             value_oppose(*d, *d);
 156             value_set_si(M->p[2*i][0], -1);
 157             Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
 158             value_set_si(M->p[2*i][0], 1);
 159             value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
 160             value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
 161         }
 162     }
 163     return M;
 164 }
 165
 166 static bool evalue_is_periodic(const evalue *e, Vector *periods)
 167 {
 168     int i, offset;
 169     bool is_periodic = false;
 170
 171     if (value_notzero_p(e->d))
 172         return false;
 173
 174     assert(e->x.p->type != partition);
 175     if (e->x.p->type == periodic) {
 176         Value size;
 177         value_init(size);
 178         value_set_si(size, e->x.p->size);
 179         value_lcm(periods->p[e->x.p->pos-1], periods->p[e->x.p->pos-1], size);
 180         value_clear(size);
 181         is_periodic = true;
 182     }
 183     offset = type_offset(e->x.p);
 184     for (i = e->x.p->size-1; i >= offset; --i)
 185         is_periodic = evalue_is_periodic(&e->x.p->arr[i], periods) || is_periodic;
 186     return is_periodic;
 187 }
 188
 189 static ex evalue2lst(const evalue *e, const exvector& vars,
 190                      exvector& extravar, vector<typed_evalue>& expr,
 191                      Vector *periods)
 192 {
 193     Vector *coset = Vector_Alloc(periods->Size);
 194     lst list;
 195     for (;;) {
 196         int i;
 197         list.append(evalue2ex_r(e, vars, extravar, expr, coset));
 198         for (i = coset->Size-1; i >= 0; --i) {
 199             value_increment(coset->p[i], coset->p[i]);
 200             if (value_lt(coset->p[i], periods->p[i]))
 201                 break;
 202             value_set_si(coset->p[i], 0);
 203         }
 204         if (i < 0)
 205             break;
 206     }
 207     Vector_Free(coset);
 208     return list;
 209 }
 210
 211 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar,
 212              Matrix **C, Vector **p)
 213 {
 214     vector<typed_evalue> expr;
 215     Vector *periods = Vector_Alloc(vars.size());
 216     assert(p);
 217     assert(C);
 218     for (int i = 0; i < periods->Size; ++i)
 219         value_set_si(periods->p[i], 1);
 220     if (evalue_is_periodic(e, periods)) {
 221         *p = periods;
 222         *C = NULL;
 223         lst list;
 224         return list;
 225     } else {
 226         Vector_Free(periods);
 227         *p = NULL;
 228         ex poly = evalue2ex_r(e, vars, floorvar, expr, NULL);
 229         Matrix *M = setup_constraints(expr, vars.size());
 230         *C = M;
 231         return poly;
 232     }
 233 }
 234
 235 /* if the evalue is a relation, we use the relation to cut off the
 236  * the edges of the domain
 237  */
 238 static Polyhedron *relation_domain(Polyhedron *D, evalue *fr, unsigned MaxRays)
 239 {
 240     assert(value_zero_p(fr->d));
 241     assert(fr->x.p->type == fractional);
 242     assert(fr->x.p->size == 3);
 243     Matrix *T = Matrix_Alloc(2, D->Dimension+1);
 244     value_set_si(T->p[1][D->Dimension], 1);
 245
 246     /* convert argument of fractional to polylib */
 247     /* the argument is assumed to be linear */
 248     evalue *p = &fr->x.p->arr[0];
 249     evalue_denom(p, &T->p[1][D->Dimension]);
 250     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
 251         assert(p->x.p->type == polynomial);
 252         assert(p->x.p->size == 2);
 253         assert(value_notzero_p(p->x.p->arr[1].d));
 254         int pos = p->x.p->pos - 1;
 255         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
 256         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][D->Dimension]);
 257         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
 258     }
 259     int pos = D->Dimension;
 260     value_assign(T->p[0][pos], p->x.n);
 261     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][D->Dimension]);
 262     value_division(T->p[0][pos], T->p[0][pos], p->d);
 263
 264     Polyhedron *E = NULL;
 265     for (Polyhedron *P = D; P; P = P->next) {
 266         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
 267         I = DomainConstraintSimplify(I, MaxRays);
 268         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
 269         Polyhedron_Free(I);
 270         Polyhedron *next = P->next;
 271         P->next = NULL;
 272         Polyhedron *S = DomainIntersection(P, R, MaxRays);
 273         Polyhedron_Free(R);
 274         P->next = next;
 275         if (emptyQ2(S))
 276             Polyhedron_Free(S);
 277         else
 278             E = DomainConcat(S, E);
 279     }
 280     Matrix_Free(T);
 281
 282     return E;
 283 }
 284
 285 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 286                                       Polyhedron *ctx, const exvector& params)
 287 {
 288     piecewise_lst *pl;
 289     barvinok_options *options = barvinok_options_new_with_defaults();
 290     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
 291     barvinok_options_free(options);
 292     return pl;
 293 }
 294
 295 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 296                             Polyhedron *D, const ex& poly,
 297                             Polyhedron *ctx,
 298                             const exvector& params, const exvector& floorvar,
 299                             barvinok_options *options);
 300
 301 /* Recursively apply Bernstein expansion on P, optimizing over dims[i]
 302  * variables in each level.  The context ctx is assumed to have been adapted
 303  * to the first level in the recursion.
 304  */
 305 static piecewise_lst *bernstein_coefficients_recursive(piecewise_lst *pl_all,
 306                             Polyhedron *P, const vector<int>& dims, const ex& poly,
 307                             Polyhedron *ctx,
 308                             const exvector& params, const exvector& vars,
 309                             barvinok_options *options)
 310 {
 311     assert(dims.size() > 0);
 312     assert(ctx->Dimension == P->Dimension - dims[0]);
 313     piecewise_lst *pl;
 314     unsigned done = 0;
 315     for (int j = 0; j < dims.size(); ++j) {
 316         exvector pl_vars;
 317         pl_vars.insert(pl_vars.end(), vars.begin()+done, vars.begin()+done+dims[j]);
 318         exvector pl_params;
 319         pl_params.insert(pl_params.end(), vars.begin()+done+dims[j], vars.end());
 320         pl_params.insert(pl_params.end(), params.begin(), params.end());
 321
 322         if (!j)
 323             pl = bernstein_coefficients(NULL, P, poly, ctx,
 324                                         pl_params, pl_vars, options);
 325         else {
 326             piecewise_lst *new_pl = NULL;
 327             Polyhedron *U = Universe_Polyhedron(pl_params.size());
 328
 329             for (int i = 0; i < pl->list.size(); ++i) {
 330                 Polyhedron *D = pl->list[i].first;
 331                 lst polys = pl->list[i].second;
 332                 new_pl = bernstein_coefficients(new_pl, D, polys, U, pl_params,
 333                                                 pl_vars, options);
 334             }
 335
 336             Polyhedron_Free(U);
 337
 338             delete pl;
 339             pl = new_pl;
 340         }
 341
 342         done += dims[j];
 343     }
 344
 345     if (!pl_all)
 346         pl_all = pl;
 347     else {
 348         pl_all->combine(*pl);
 349         delete pl;
 350     }
 351
 352     return pl_all;
 353 }
 354
 355 static piecewise_lst *bernstein_coefficients_full_recurse(piecewise_lst *pl_all,
 356                             Polyhedron *P, const ex& poly,
 357                             Polyhedron *ctx,
 358                             const exvector& params, const exvector& vars,
 359                             barvinok_options *options)
 360 {
 361     Polyhedron *CR = align_context(ctx, P->Dimension-1, options->MaxRays);
 362     vector<int> dims(vars.size());
 363     for (int i = 0; i < dims.size(); ++i)
 364         dims[i] = 1;
 365     pl_all = bernstein_coefficients_recursive(pl_all, P, dims, poly, CR,
 366                                               params, vars, options);
 367     Polyhedron_Free(CR);
 368
 369     return pl_all;
 370 }
 371
 372 static piecewise_lst *bernstein_coefficients_product(piecewise_lst *pl_all,
 373                             Polyhedron *F, Matrix *T, const ex& poly,
 374                             Polyhedron *ctx,
 375                             const exvector& params, const exvector& vars,
 376                             barvinok_options *options)
 377 {
 378     if (emptyQ2(ctx))
 379         return pl_all;
 380     for (Polyhedron *G = F; G; G = G->next)
 381         if (emptyQ2(G))
 382             return pl_all;
 383
 384     unsigned nparam = params.size();
 385     unsigned nvar = vars.size();
 386     unsigned constraints;
 387     unsigned factors;
 388     Polyhedron *C = NULL;
 389
 390     /* More context constraints */
 391     if (F->Dimension == ctx->Dimension) {
 392         C = F;
 393         F = F->next;
 394     }
 395     assert(F);
 396     assert(F->next);
 397
 398     Matrix *M;
 399     Polyhedron *P;
 400     Polyhedron *PC;
 401     M = Matrix_Alloc(F->NbConstraints, 1+nvar+nparam+1);
 402     for (int i = 0; i < F->NbConstraints; ++i) {
 403         Vector_Copy(F->Constraint[i], M->p[i], 1+F->Dimension-nparam);
 404         Vector_Copy(F->Constraint[i]+1+F->Dimension-nparam,
 405                     M->p[i]+1+nvar, nparam+1);
 406     }
 407     P = Constraints2Polyhedron(M, options->MaxRays);
 408     Matrix_Free(M);
 409
 410     factors = 1;
 411     constraints = C ? C->NbConstraints : 0;
 412     constraints += ctx->NbConstraints;
 413     for (Polyhedron *G = F->next; G; G = G->next) {
 414         constraints += G->NbConstraints;
 415         ++factors;
 416     }
 417
 418     unsigned total_var = nvar-(F->Dimension-nparam);
 419     unsigned skip = 0;
 420     unsigned c = 0;
 421     M = Matrix_Alloc(constraints, 1+total_var+nparam+1);
 422     for (Polyhedron *G = F->next; G; G = G->next) {
 423         unsigned this_var = G->Dimension - nparam;
 424         for (int i = 0; i < G->NbConstraints; ++i) {
 425             value_assign(M->p[c+i][0], G->Constraint[i][0]);
 426             Vector_Copy(G->Constraint[i]+1, M->p[c+i]+1+skip, this_var);
 427             Vector_Copy(G->Constraint[i]+1+this_var, M->p[c+i]+1+total_var,
 428                         nparam+1);
 429         }
 430         c += G->NbConstraints;
 431         skip += this_var;
 432     }
 433     assert(skip == total_var);
 434     if (C) {
 435         for (int i = 0; i < C->NbConstraints; ++i) {
 436             value_assign(M->p[c+i][0], C->Constraint[i][0]);
 437             Vector_Copy(C->Constraint[i]+1, M->p[c+i]+1+total_var,
 438                         nparam+1);
 439         }
 440         c += C->NbConstraints;
 441     }
 442     for (int i = 0; i < ctx->NbConstraints; ++i) {
 443         value_assign(M->p[c+i][0], ctx->Constraint[i][0]);
 444         Vector_Copy(ctx->Constraint[i]+1, M->p[c+i]+1+total_var, nparam+1);
 445     }
 446     PC = Constraints2Polyhedron(M, options->MaxRays);
 447     Matrix_Free(M);
 448
 449     exvector newvars = constructVariableVector(nvar, "t");
 450     matrix subs(1, nvar);
 451     for (int i = 0; i < nvar; ++i)
 452         for (int j = 0; j < nvar; ++j)
 453             subs(0,i) += value2numeric(T->p[i][j]) * newvars[j];
 454
 455     ex newpoly = replaceVariablesInPolynomial(poly, vars, subs);
 456
 457     vector<int> dims(factors);
 458     for (int i = 0; F; ++i, F = F->next)
 459         dims[i] = F->Dimension-nparam;
 460
 461     pl_all = bernstein_coefficients_recursive(pl_all, P, dims, newpoly, PC,
 462                                               params, newvars, options);
 463
 464     Polyhedron_Free(P);
 465     Polyhedron_Free(PC);
 466
 467     return pl_all;
 468 }
 469
 470 static piecewise_lst *bernstein_coefficients_polyhedron(piecewise_lst *pl_all,
 471                             Polyhedron *P, const ex& poly,
 472                             Polyhedron *ctx,
 473                             const exvector& params, const exvector& floorvar,
 474                             barvinok_options *options)
 475 {
 476     if (Polyhedron_is_unbounded(P, ctx->Dimension, options->MaxRays)) {
 477         fprintf(stderr, "warning: unbounded domain skipped\n");
 478         Polyhedron_Print(stderr, P_VALUE_FMT, P);
 479         return pl_all;
 480     }
 481
 482     if (options->bernstein_recurse & BV_BERNSTEIN_FACTORS) {
 483         Matrix *T = NULL;
 484         Polyhedron *F = Polyhedron_Factor(P, ctx->Dimension, &T, options->MaxRays);
 485         if (F) {
 486             pl_all = bernstein_coefficients_product(pl_all, F, T, poly, ctx, params,
 487                                                     floorvar, options);
 488             Domain_Free(F);
 489             Matrix_Free(T);
 490             return pl_all;
 491         }
 492     }
 493     if (floorvar.size() > 1 &&
 494                 options->bernstein_recurse & BV_BERNSTEIN_INTERVALS)
 495         return bernstein_coefficients_full_recurse(pl_all, P, poly, ctx,
 496                                                    params, floorvar, options);
 497
 498     unsigned PP_MaxRays = options->MaxRays;
 499     if (PP_MaxRays & POL_NO_DUAL)
 500         PP_MaxRays = 0;
 501
 502     Param_Polyhedron *PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 503     assert(PP);
 504     piecewise_lst *pl = new piecewise_lst(params, options->bernstein_optimize);
 505     for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 506         matrix VM = domainVertices(PP, Q, params);
 507         lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 508         pl->add_guarded_lst(Polyhedron_Copy(Q->Domain), coeffs);
 509     }
 510     Param_Polyhedron_Free(PP);
 511     if (!pl_all)
 512         pl_all = pl;
 513     else {
 514         pl_all->combine(*pl);
 515         delete pl;
 516     }
 517
 518     return pl_all;
 519 }
 520
 521 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 522                             Polyhedron *D, const ex& poly,
 523                             Polyhedron *ctx,
 524                             const exvector& params, const exvector& floorvar,
 525                             barvinok_options *options)
 526 {
 527     if (!D->next && emptyQ2(D))
 528         return pl_all;
 529
 530     for (Polyhedron *P = D; P; P = P->next) {
 531         /* This shouldn't happen */
 532         if (emptyQ2(P))
 533             continue;
 534         Polyhedron *next = P->next;
 535         P->next = NULL;
 536         pl_all = bernstein_coefficients_polyhedron(pl_all, P, poly, ctx,
 537                                                    params, floorvar, options);
 538         P->next = next;
 539     }
 540     return pl_all;
 541 }
 542
 543 /* Compute the coefficients of the polynomial corresponding to each coset
 544  * on its own domain.  This allows us to cut the domain on multiples of
 545  * the period.
 546  * To perform the cutting for a coset "i mod n = c" we map the domain
 547  * to the quotient space trough "i = i' n + c", simplify the constraints
 548  * (implicitly) and then map back to the original space.
 549  */
 550 static piecewise_lst *bernstein_coefficients_periodic(piecewise_lst *pl_all,
 551                             Polyhedron *D, const evalue *e,
 552                             Polyhedron *ctx, const exvector& vars,
 553                             const exvector& params, Vector *periods,
 554                             barvinok_options *options)
 555 {
 556     assert(D->Dimension == periods->Size);
 557     Matrix *T = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 558     Matrix *T2 = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 559     Vector *coset = Vector_Alloc(periods->Size);
 560     exvector extravar;
 561     vector<typed_evalue> expr;
 562     exvector allvars = vars;
 563     allvars.insert(allvars.end(), params.begin(), params.end());
 564
 565     value_set_si(T2->p[D->Dimension][D->Dimension], 1);
 566     for (int i = 0; i < D->Dimension; ++i) {
 567         value_assign(T->p[i][i], periods->p[i]);
 568         value_lcm(T2->p[D->Dimension][D->Dimension],
 569                   T2->p[D->Dimension][D->Dimension], periods->p[i]);
 570     }
 571     value_set_si(T->p[D->Dimension][D->Dimension], 1);
 572     for (int i = 0; i < D->Dimension; ++i)
 573         value_division(T2->p[i][i], T2->p[D->Dimension][D->Dimension],
 574                        periods->p[i]);
 575     for (;;) {
 576         int i;
 577         ex poly = evalue2ex_r(e, allvars, extravar, expr, coset);
 578         assert(extravar.size() == 0);
 579         assert(expr.size() == 0);
 580         Polyhedron *E = DomainPreimage(D, T, options->MaxRays);
 581         Polyhedron *F = DomainPreimage(E, T2, options->MaxRays);
 582         Polyhedron_Free(E);
 583         if (!emptyQ2(F))
 584             pl_all = bernstein_coefficients(pl_all, F, poly, ctx, params,
 585                                             vars, options);
 586         Polyhedron_Free(F);
 587         for (i = D->Dimension-1; i >= 0; --i) {
 588             value_increment(coset->p[i], coset->p[i]);
 589             value_increment(T->p[i][D->Dimension], T->p[i][D->Dimension]);
 590             value_subtract(T2->p[i][D->Dimension], T2->p[i][D->Dimension],
 591                            T2->p[i][i]);
 592             if (value_lt(coset->p[i], periods->p[i]))
 593                 break;
 594             value_set_si(coset->p[i], 0);
 595             value_set_si(T->p[i][D->Dimension], 0);
 596             value_set_si(T2->p[i][D->Dimension], 0);
 597         }
 598         if (i < 0)
 599             break;
 600     }
 601     Vector_Free(coset);
 602     Matrix_Free(T);
 603     Matrix_Free(T2);
 604     return pl_all;
 605 }
 606
 607 piecewise_lst *bernstein_coefficients_relation(piecewise_lst *pl_all,
 608                     Polyhedron *D, evalue *EP, Polyhedron *ctx,
 609                     const exvector& allvars, const exvector& vars,
 610                     const exvector& params, barvinok_options *options)
 611 {
 612     if (value_zero_p(EP->d) && EP->x.p->type == relation) {
 613         Polyhedron *E = relation_domain(D, &EP->x.p->arr[0], options->MaxRays);
 614         if (E) {
 615             pl_all = bernstein_coefficients_relation(pl_all, E, &EP->x.p->arr[1],
 616                                                      ctx, allvars, vars, params,
 617                                                      options);
 618             Domain_Free(E);
 619         }
 620         /* In principle, we could cut off the edges of this domain too */
 621         if (EP->x.p->size > 2)
 622             pl_all = bernstein_coefficients_relation(pl_all, D, &EP->x.p->arr[2],
 623                                                      ctx, allvars, vars, params,
 624                                                      options);
 625         return pl_all;
 626     }
 627
 628     Matrix *M;
 629     exvector floorvar;
 630     Vector *periods;
 631     ex poly = evalue2ex(EP, allvars, floorvar, &M, &periods);
 632     floorvar.insert(floorvar.end(), vars.begin(), vars.end());
 633     Polyhedron *E = D;
 634     if (M) {
 635         Polyhedron *AE = align_context(D, M->NbColumns-2, options->MaxRays);
 636         E = DomainAddConstraints(AE, M, options->MaxRays);
 637         Matrix_Free(M);
 638         Domain_Free(AE);
 639     }
 640     if (is_exactly_a<fail>(poly)) {
 641         delete pl_all;
 642         return NULL;
 643     }
 644     if (periods)
 645         pl_all = bernstein_coefficients_periodic(pl_all, E, EP, ctx, vars,
 646                                                  params, periods, options);
 647     else
 648         pl_all = bernstein_coefficients(pl_all, E, poly, ctx, params,
 649                                         floorvar, options);
 650     if (periods)
 651         Vector_Free(periods);
 652     if (D != E)
 653         Domain_Free(E);
 654
 655     return pl_all;
 656 }
 657
 658 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 659                                       Polyhedron *ctx, const exvector& params,
 660                                       barvinok_options *options)
 661 {
 662     unsigned nparam = ctx->Dimension;
 663     if (EVALUE_IS_ZERO(*e))
 664         return pl_all;
 665     assert(value_zero_p(e->d));
 666     assert(e->x.p->type == partition);
 667     assert(e->x.p->size >= 2);
 668     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 669
 670     exvector vars = constructVariableVector(nvars, "v");
 671     exvector allvars = vars;
 672     allvars.insert(allvars.end(), params.begin(), params.end());
 673
 674     for (int i = 0; i < e->x.p->size/2; ++i) {
 675         pl_all = bernstein_coefficients_relation(pl_all,
 676                         EVALUE_DOMAIN(e->x.p->arr[2*i]), &e->x.p->arr[2*i+1],
 677                         ctx, allvars, vars, params, options);
 678     }
 679     return pl_all;
 680 }
 681
 682 }