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             value_set_si(*d, 1);
 153             evalue_denom(expr[i].second, d);
 154             const evalue *e;
 155             for (e = expr[i].second; value_zero_p(e->d); e = &e->x.p->arr[0]) {
 156                 assert(e->x.p->type == polynomial);
 157                 assert(e->x.p->size == 2);
 158                 evalue *c = &e->x.p->arr[1];
 159                 value_multiply(M->p[2*i][1+extra+e->x.p->pos-1], *d, c->x.n);
 160                 value_division(M->p[2*i][1+extra+e->x.p->pos-1],
 161                                M->p[2*i][1+extra+e->x.p->pos-1], c->d);
 162             }
 163             value_multiply(M->p[2*i][1+extra+nvar], *d, e->x.n);
 164             value_division(M->p[2*i][1+extra+nvar], M->p[2*i][1+extra+nvar], e->d);
 165             value_oppose(*d, *d);
 166             value_set_si(M->p[2*i][0], -1);
 167             Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
 168             value_set_si(M->p[2*i][0], 1);
 169             value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
 170             value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
 171         }
 172     }
 173     return M;
 174 }
 175
 176 static bool evalue_is_periodic(const evalue *e, Vector *periods)
 177 {
 178     int i, offset;
 179     bool is_periodic = false;
 180
 181     if (value_notzero_p(e->d))
 182         return false;
 183
 184     assert(e->x.p->type != partition);
 185     if (e->x.p->type == periodic) {
 186         Value size;
 187         value_init(size);
 188         value_set_si(size, e->x.p->size);
 189         value_lcm(periods->p[e->x.p->pos-1], size, &periods->p[e->x.p->pos-1]);
 190         value_clear(size);
 191         is_periodic = true;
 192     }
 193     offset = type_offset(e->x.p);
 194     for (i = e->x.p->size-1; i >= offset; --i)
 195         is_periodic = evalue_is_periodic(&e->x.p->arr[i], periods) || is_periodic;
 196     return is_periodic;
 197 }
 198
 199 static ex evalue2lst(const evalue *e, const exvector& vars,
 200                      exvector& extravar, vector<typed_evalue>& expr,
 201                      Vector *periods)
 202 {
 203     Vector *coset = Vector_Alloc(periods->Size);
 204     lst list;
 205     for (;;) {
 206         int i;
 207         list.append(evalue2ex_r(e, vars, extravar, expr, coset));
 208         for (i = coset->Size-1; i >= 0; --i) {
 209             value_increment(coset->p[i], coset->p[i]);
 210             if (value_lt(coset->p[i], periods->p[i]))
 211                 break;
 212             value_set_si(coset->p[i], 0);
 213         }
 214         if (i < 0)
 215             break;
 216     }
 217     Vector_Free(coset);
 218     return list;
 219 }
 220
 221 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar,
 222              Matrix **C, Vector **p)
 223 {
 224     vector<typed_evalue> expr;
 225     Vector *periods = Vector_Alloc(vars.size());
 226     assert(p);
 227     assert(C);
 228     for (int i = 0; i < periods->Size; ++i)
 229         value_set_si(periods->p[i], 1);
 230     if (evalue_is_periodic(e, periods)) {
 231         *p = periods;
 232         *C = NULL;
 233         lst list;
 234         return list;
 235     } else {
 236         Vector_Free(periods);
 237         *p = NULL;
 238         ex poly = evalue2ex_r(e, vars, floorvar, expr, NULL);
 239         Matrix *M = setup_constraints(expr, vars.size());
 240         *C = M;
 241         return poly;
 242     }
 243 }
 244
 245 /* if the evalue is a relation, we use the relation to cut off the
 246  * the edges of the domain
 247  */
 248 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
 249                                 unsigned MaxRays)
 250 {
 251     *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 252     *poly = e = &e->x.p->arr[2*i+1];
 253     if (value_notzero_p(e->d))
 254         return;
 255     if (e->x.p->type != relation)
 256         return;
 257     if (e->x.p->size > 2)
 258         return;
 259     evalue *fr = &e->x.p->arr[0];
 260     assert(value_zero_p(fr->d));
 261     assert(fr->x.p->type == fractional);
 262     assert(fr->x.p->size == 3);
 263     Matrix *T = Matrix_Alloc(2, (*D)->Dimension+1);
 264     value_set_si(T->p[1][(*D)->Dimension], 1);
 265
 266     /* convert argument of fractional to polylib */
 267     /* the argument is assumed to be linear */
 268     evalue *p = &fr->x.p->arr[0];
 269     evalue_denom(p, &T->p[1][(*D)->Dimension]);
 270     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
 271         assert(p->x.p->type == polynomial);
 272         assert(p->x.p->size == 2);
 273         assert(value_notzero_p(p->x.p->arr[1].d));
 274         int pos = p->x.p->pos - 1;
 275         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
 276         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 277         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
 278     }
 279     int pos = (*D)->Dimension;
 280     value_assign(T->p[0][pos], p->x.n);
 281     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 282     value_division(T->p[0][pos], T->p[0][pos], p->d);
 283
 284     Polyhedron *E = NULL;
 285     for (Polyhedron *P = *D; P; P = P->next) {
 286         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
 287         I = DomainConstraintSimplify(I, MaxRays);
 288         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
 289         Polyhedron_Free(I);
 290         Polyhedron *next = P->next;
 291         P->next = NULL;
 292         Polyhedron *S = DomainIntersection(P, R, MaxRays);
 293         Polyhedron_Free(R);
 294         P->next = next;
 295         if (emptyQ2(S))
 296             Polyhedron_Free(S);
 297         else
 298             E = DomainConcat(S, E);
 299     }
 300     Matrix_Free(T);
 301
 302     *D = E;
 303     *poly = &e->x.p->arr[1];
 304 }
 305
 306 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 307                                       Polyhedron *ctx, const exvector& params)
 308 {
 309     piecewise_lst *pl;
 310     barvinok_options *options = barvinok_options_new_with_defaults();
 311     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
 312     barvinok_options_free(options);
 313     return pl;
 314 }
 315
 316 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 317                             Polyhedron *D, const ex& poly,
 318                             Polyhedron *ctx,
 319                             const exvector& params, const exvector& floorvar,
 320                             barvinok_options *options);
 321
 322 static piecewise_lst *bernstein_coefficients_product(piecewise_lst *pl_all,
 323                             Polyhedron *F, Matrix *T, const ex& poly,
 324                             Polyhedron *ctx,
 325                             const exvector& params, const exvector& vars,
 326                             barvinok_options *options)
 327 {
 328     if (emptyQ2(ctx))
 329         return pl_all;
 330     for (Polyhedron *G = F; G; G = G->next)
 331         if (emptyQ2(G))
 332             return pl_all;
 333
 334     unsigned nparam = params.size();
 335     unsigned nvar = vars.size();
 336     unsigned constraints;
 337     Polyhedron *C = NULL;
 338
 339     /* More context constraints */
 340     if (F->Dimension == ctx->Dimension) {
 341         C = F;
 342         F = F->next;
 343     }
 344     assert(F);
 345     assert(F->next);
 346
 347     Matrix *M;
 348     Polyhedron *P;
 349     Polyhedron *PC;
 350     M = Matrix_Alloc(F->NbConstraints, 1+nvar+nparam+1);
 351     for (int i = 0; i < F->NbConstraints; ++i) {
 352         Vector_Copy(F->Constraint[i], M->p[i], 1+F->Dimension-nparam);
 353         Vector_Copy(F->Constraint[i]+1+F->Dimension-nparam,
 354                     M->p[i]+1+nvar, nparam+1);
 355     }
 356     P = Constraints2Polyhedron(M, options->MaxRays);
 357     Matrix_Free(M);
 358
 359     constraints = C ? C->NbConstraints : 0;
 360     constraints += ctx->NbConstraints;
 361     for (Polyhedron *G = F->next; G; G = G->next)
 362         constraints += G->NbConstraints;
 363
 364     unsigned total_var = nvar-(F->Dimension-nparam);
 365     unsigned skip = 0;
 366     unsigned c = 0;
 367     M = Matrix_Alloc(constraints, 1+total_var+nparam+1);
 368     for (Polyhedron *G = F->next; G; G = G->next) {
 369         unsigned this_var = G->Dimension - nparam;
 370         for (int i = 0; i < G->NbConstraints; ++i) {
 371             value_assign(M->p[c+i][0], G->Constraint[i][0]);
 372             Vector_Copy(G->Constraint[i]+1, M->p[c+i]+1+skip, this_var);
 373             Vector_Copy(G->Constraint[i]+1+this_var, M->p[c+i]+1+total_var,
 374                         nparam+1);
 375         }
 376         c += G->NbConstraints;
 377         skip += this_var;
 378     }
 379     assert(skip == total_var);
 380     if (C) {
 381         for (int i = 0; i < C->NbConstraints; ++i) {
 382             value_assign(M->p[c+i][0], C->Constraint[i][0]);
 383             Vector_Copy(C->Constraint[i]+1, M->p[c+i]+1+total_var,
 384                         nparam+1);
 385         }
 386         c += C->NbConstraints;
 387     }
 388     for (int i = 0; i < ctx->NbConstraints; ++i) {
 389         value_assign(M->p[c+i][0], ctx->Constraint[i][0]);
 390         Vector_Copy(ctx->Constraint[i]+1, M->p[c+i]+1+total_var, nparam+1);
 391     }
 392     PC = Constraints2Polyhedron(M, options->MaxRays);
 393     Matrix_Free(M);
 394
 395     exvector newvars = constructVariableVector(nvar, "t");
 396     matrix subs(1, nvar);
 397     for (int i = 0; i < nvar; ++i)
 398         for (int j = 0; j < nvar; ++j)
 399             subs(0,i) += value2numeric(T->p[i][j]) * newvars[j];
 400
 401     ex newpoly = replaceVariablesInPolynomial(poly, vars, subs);
 402
 403     exvector P_vars;
 404     P_vars.insert(P_vars.end(), newvars.begin(),
 405                   newvars.begin()+(F->Dimension-nparam));
 406     exvector P_params;
 407     P_params.insert(P_params.end(), newvars.begin()+(F->Dimension-nparam),
 408                     newvars.end());
 409     P_params.insert(P_params.end(), params.begin(), params.end());
 410     piecewise_lst *pl;
 411     pl = bernstein_coefficients(NULL, P, newpoly, PC, P_params, P_vars, options);
 412     Polyhedron_Free(P);
 413     Polyhedron_Free(PC);
 414
 415     unsigned done = F->Dimension-nparam;
 416     for (F = F->next ; F; F = F->next) {
 417         exvector pl_vars;
 418         pl_vars.insert(pl_vars.end(), newvars.begin()+done,
 419                        newvars.begin()+done+(F->Dimension-nparam));
 420         exvector pl_params;
 421         pl_params.insert(pl_params.end(), newvars.begin()+done+(F->Dimension-nparam),
 422                          newvars.end());
 423         pl_params.insert(pl_params.end(), params.begin(), params.end());
 424         piecewise_lst *new_pl = NULL;
 425         Polyhedron *U = Universe_Polyhedron(pl_params.size());
 426
 427         for (int i = 0; i < pl->list.size(); ++i) {
 428             Polyhedron *D = pl->list[i].first;
 429             lst polys = pl->list[i].second;
 430             lst::const_iterator j;
 431             for (j = polys.begin(); j != polys.end(); ++j) {
 432                 new_pl = bernstein_coefficients(new_pl, D, *j, U, pl_params,
 433                                                 pl_vars, options);
 434             }
 435         }
 436
 437         delete pl;
 438         pl = new_pl;
 439
 440         Polyhedron_Free(U);
 441         done += F->Dimension-nparam;
 442     }
 443
 444     if (!pl_all)
 445         pl_all = pl;
 446     else {
 447         pl_all->combine(*pl);
 448         delete pl;
 449     }
 450
 451     return pl_all;
 452 }
 453
 454 static piecewise_lst *bernstein_coefficients_polyhedron(piecewise_lst *pl_all,
 455                             Polyhedron *P, const ex& poly,
 456                             Polyhedron *ctx,
 457                             const exvector& params, const exvector& floorvar,
 458                             barvinok_options *options)
 459 {
 460     if (Polyhedron_is_unbounded(P, ctx->Dimension, options->MaxRays)) {
 461         fprintf(stderr, "warning: unbounded domain skipped\n");
 462         Polyhedron_Print(stderr, P_VALUE_FMT, P);
 463         return pl_all;
 464     }
 465
 466     Matrix *T = NULL;
 467     Polyhedron *F = Polyhedron_Factor(P, ctx->Dimension, &T, options->MaxRays);
 468     if (F) {
 469         pl_all = bernstein_coefficients_product(pl_all, F, T, poly, ctx, params,
 470                                                 floorvar, options);
 471         Domain_Free(F);
 472         Matrix_Free(T);
 473         return pl_all;
 474     }
 475
 476     unsigned PP_MaxRays = options->MaxRays;
 477     if (PP_MaxRays & POL_NO_DUAL)
 478         PP_MaxRays = 0;
 479
 480     Param_Polyhedron *PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 481     assert(PP);
 482     piecewise_lst *pl = new piecewise_lst(params);
 483     for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 484         matrix VM = domainVertices(PP, Q, params);
 485         lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 486         pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
 487     }
 488     Param_Polyhedron_Free(PP);
 489     if (!pl_all)
 490         pl_all = pl;
 491     else {
 492         pl_all->combine(*pl);
 493         delete pl;
 494     }
 495
 496     return pl_all;
 497 }
 498
 499 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 500                             Polyhedron *D, const ex& poly,
 501                             Polyhedron *ctx,
 502                             const exvector& params, const exvector& floorvar,
 503                             barvinok_options *options)
 504 {
 505     if (!D->next && emptyQ2(D))
 506         return pl_all;
 507
 508     for (Polyhedron *P = D; P; P = P->next) {
 509         /* This shouldn't happen */
 510         if (emptyQ2(P))
 511             continue;
 512         Polyhedron *next = P->next;
 513         P->next = NULL;
 514         pl_all = bernstein_coefficients_polyhedron(pl_all, P, poly, ctx,
 515                                                    params, floorvar, options);
 516         P->next = next;
 517     }
 518     return pl_all;
 519 }
 520
 521 /* Compute the coefficients of the polynomial corresponding to each coset
 522  * on its own domain.  This allows us to cut the domain on multiples of
 523  * the period.
 524  * To perform the cutting for a coset "i mod n = c" we map the domain
 525  * to the quotient space trough "i = i' n + c", simplify the constraints
 526  * (implicitly) and then map back to the original space.
 527  */
 528 static piecewise_lst *bernstein_coefficients_periodic(piecewise_lst *pl_all,
 529                             Polyhedron *D, const evalue *e,
 530                             Polyhedron *ctx, const exvector& vars,
 531                             const exvector& params, Vector *periods,
 532                             barvinok_options *options)
 533 {
 534     assert(D->Dimension == periods->Size);
 535     Matrix *T = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 536     Matrix *T2 = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 537     Vector *coset = Vector_Alloc(periods->Size);
 538     exvector extravar;
 539     vector<typed_evalue> expr;
 540     exvector allvars = vars;
 541     allvars.insert(allvars.end(), params.begin(), params.end());
 542
 543     value_set_si(T2->p[D->Dimension][D->Dimension], 1);
 544     for (int i = 0; i < D->Dimension; ++i) {
 545         value_assign(T->p[i][i], periods->p[i]);
 546         value_lcm(T2->p[D->Dimension][D->Dimension], periods->p[i],
 547                   &T2->p[D->Dimension][D->Dimension]);
 548     }
 549     value_set_si(T->p[D->Dimension][D->Dimension], 1);
 550     for (int i = 0; i < D->Dimension; ++i)
 551         value_division(T2->p[i][i], T2->p[D->Dimension][D->Dimension],
 552                        periods->p[i]);
 553     for (;;) {
 554         int i;
 555         ex poly = evalue2ex_r(e, allvars, extravar, expr, coset);
 556         assert(extravar.size() == 0);
 557         assert(expr.size() == 0);
 558         Polyhedron *E = DomainPreimage(D, T, options->MaxRays);
 559         Polyhedron *F = DomainPreimage(E, T2, options->MaxRays);
 560         Polyhedron_Free(E);
 561         if (!emptyQ2(F))
 562             pl_all = bernstein_coefficients(pl_all, F, poly, ctx, params,
 563                                             vars, options);
 564         Polyhedron_Free(F);
 565         for (i = D->Dimension-1; i >= 0; --i) {
 566             value_increment(coset->p[i], coset->p[i]);
 567             value_increment(T->p[i][D->Dimension], T->p[i][D->Dimension]);
 568             value_subtract(T2->p[i][D->Dimension], T2->p[i][D->Dimension],
 569                            T2->p[i][i]);
 570             if (value_lt(coset->p[i], periods->p[i]))
 571                 break;
 572             value_set_si(coset->p[i], 0);
 573             value_set_si(T->p[i][D->Dimension], 0);
 574             value_set_si(T2->p[i][D->Dimension], 0);
 575         }
 576         if (i < 0)
 577             break;
 578     }
 579     Vector_Free(coset);
 580     Matrix_Free(T);
 581     Matrix_Free(T2);
 582     return pl_all;
 583 }
 584
 585 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 586                                       Polyhedron *ctx, const exvector& params,
 587                                       barvinok_options *options)
 588 {
 589     unsigned nparam = ctx->Dimension;
 590     if (EVALUE_IS_ZERO(*e))
 591         return pl_all;
 592     assert(value_zero_p(e->d));
 593     assert(e->x.p->type == partition);
 594     assert(e->x.p->size >= 2);
 595     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 596
 597     exvector vars = constructVariableVector(nvars, "v");
 598     exvector allvars = vars;
 599     allvars.insert(allvars.end(), params.begin(), params.end());
 600
 601     for (int i = 0; i < e->x.p->size/2; ++i) {
 602         Polyhedron *E;
 603         evalue *EP;
 604         Matrix *M;
 605         Vector *periods;
 606         exvector floorvar;
 607
 608         evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
 609         ex poly = evalue2ex(EP, allvars, floorvar, &M, &periods);
 610         floorvar.insert(floorvar.end(), vars.begin(), vars.end());
 611         if (M) {
 612             Polyhedron *AE = align_context(E, M->NbColumns-2, options->MaxRays);
 613             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 614                 Domain_Free(E);
 615             E = DomainAddConstraints(AE, M, options->MaxRays);
 616             Matrix_Free(M);
 617             Domain_Free(AE);
 618         }
 619         if (is_exactly_a<fail>(poly)) {
 620             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 621                 Domain_Free(E);
 622             delete pl_all;
 623             return NULL;
 624         }
 625         if (periods)
 626             pl_all = bernstein_coefficients_periodic(pl_all, E, EP, ctx, vars,
 627                                                      params, periods, options);
 628         else
 629             pl_all = bernstein_coefficients(pl_all, E, poly, ctx, params,
 630                                             floorvar, options);
 631         if (periods)
 632             Vector_Free(periods);
 633         if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 634             Domain_Free(E);
 635     }
 636     return pl_all;
 637 }
 638
 639 }