bernstein.cc

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