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