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