bernstein.cc

   1 #include <assert.h>
   2 #include <vector>
   3 #include <bernstein/bernstein.h>
   4 #include <bernstein/piecewise_lst.h>
   5 #include <isl_set_polylib.h>
   6 #include <barvinok/barvinok.h>
   7 #include <barvinok/util.h>
   8 #include <barvinok/bernstein.h>
   9 #include <barvinok/options.h>
  10 #include "bound_common.h"
  11 #include "reduce_domain.h"
  12
  13 using namespace GiNaC;
  14 using namespace bernstein;
  15 using namespace barvinok;
  16
  17 using std::pair;
  18 using std::vector;
  19 using std::cerr;
  20 using std::endl;
  21
  22 namespace barvinok {
  23
  24 ex evalue2ex(evalue *e, const exvector& vars)
  25 {
  26     if (value_pos_p(e->d))
  27         return value2numeric(e->x.n)/value2numeric(e->d);
  28     if (EVALUE_IS_NAN(*e))
  29         return fail();
  30     if (e->x.p->type != polynomial)
  31         return fail();
  32     ex poly = 0;
  33     for (int i = e->x.p->size-1; i >= 0; --i) {
  34         poly *= vars[e->x.p->pos-1];
  35         ex t = evalue2ex(&e->x.p->arr[i], vars);
  36         if (is_exactly_a<fail>(t))
  37             return t;
  38         poly += t;
  39     }
  40     return poly;
  41 }
  42
  43 static int type_offset(enode *p)
  44 {
  45    return p->type == fractional ? 1 :
  46           p->type == flooring ? 1 : 0;
  47 }
  48
  49 typedef pair<bool, const evalue *> typed_evalue;
  50
  51 static ex evalue2ex_add_var(evalue *e, exvector& extravar,
  52                             vector<typed_evalue>& expr, bool is_fract)
  53 {
  54     ex base_var = 0;
  55
  56     for (int i = 0; i < expr.size(); ++i) {
  57         if (is_fract == expr[i].first && eequal(e, expr[i].second)) {
  58             base_var = extravar[i];
  59             break;
  60         }
  61     }
  62     if (base_var != 0)
  63         return base_var;
  64
  65     char name[20];
  66     snprintf(name, sizeof(name), "f%c%zd", is_fract ? 'r' : 'l', expr.size());
  67     extravar.push_back(base_var = symbol(name));
  68     expr.push_back(typed_evalue(is_fract, e));
  69
  70     return base_var;
  71 }
  72
  73 /* For the argument e=(f/d) of a fractional, return (d-1)/d times
  74  * a variable in [0,1] (see setup_constraints).
  75  */
  76 static ex evalue2ex_get_fract(evalue *e, exvector& extravar,
  77                                 vector<typed_evalue>& expr)
  78 {
  79     ex f;
  80     Value d;
  81     ex den;
  82     value_init(d);
  83     value_set_si(d, 1);
  84     evalue_denom(e, &d);
  85     den = value2numeric(d);
  86     value_clear(d);
  87     f = (den-1)/den;
  88
  89     ex base_var = evalue2ex_add_var(e, extravar, expr, true);
  90     base_var *= f;
  91     return base_var;
  92 }
  93
  94 static ex evalue2ex_r(const evalue *e, const exvector& vars,
  95                       exvector& extravar, vector<typed_evalue>& expr,
  96                       Vector *coset)
  97 {
  98     if (value_notzero_p(e->d))
  99         return value2numeric(e->x.n)/value2numeric(e->d);
 100     ex base_var = 0;
 101     ex poly = 0;
 102     int rem;
 103
 104     switch (e->x.p->type) {
 105     case polynomial:
 106         base_var = vars[e->x.p->pos-1];
 107         break;
 108     case flooring:
 109         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, false);
 110         break;
 111     case fractional:
 112         base_var = evalue2ex_get_fract(&e->x.p->arr[0], extravar, expr);
 113         break;
 114     case periodic:
 115         assert(coset);
 116         rem = VALUE_TO_INT(coset->p[e->x.p->pos-1]) % e->x.p->size;
 117         return evalue2ex_r(&e->x.p->arr[rem], vars, extravar, expr, coset);
 118     default:
 119         return fail();
 120     }
 121
 122     int offset = type_offset(e->x.p);
 123     for (int i = e->x.p->size-1; i >= offset; --i) {
 124         poly *= base_var;
 125         ex t = evalue2ex_r(&e->x.p->arr[i], vars, extravar, expr, coset);
 126         if (is_exactly_a<fail>(t))
 127             return t;
 128         poly += t;
 129     }
 130     return poly;
 131 }
 132
 133 /* For each t = floor(e/d), set up two constraints
 134  *
 135  *       e - d t       >= 0
 136  *      -e + d t + d-1 >= 0
 137  *
 138  * e is assumed to be an affine expression.
 139  *
 140  * For each t = fract(e/d), set up two constraints
 141  *
 142  *          -d t + d-1 >= 0
 143  *             t       >= 0
 144  */
 145 static Matrix *setup_constraints(const vector<typed_evalue> expr, int nvar)
 146 {
 147     int extra = expr.size();
 148     if (!extra)
 149         return NULL;
 150     Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
 151     for (int i = 0; i < extra; ++i) {
 152         if (expr[i].first) {
 153             value_set_si(M->p[2*i][0], 1);
 154             value_set_si(M->p[2*i][1+i], -1);
 155             value_set_si(M->p[2*i][1+extra+nvar], 1);
 156             value_set_si(M->p[2*i+1][0], 1);
 157             value_set_si(M->p[2*i+1][1+i], 1);
 158         } else {
 159             Value *d = &M->p[2*i][1+i];
 160             evalue_extract_affine(expr[i].second, M->p[2*i]+1+extra,
 161                                   M->p[2*i]+1+extra+nvar, d);
 162             value_oppose(*d, *d);
 163             value_set_si(M->p[2*i][0], -1);
 164             Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
 165             value_set_si(M->p[2*i][0], 1);
 166             value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
 167             value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
 168         }
 169     }
 170     return M;
 171 }
 172
 173 static bool evalue_is_periodic(const evalue *e, Vector *periods)
 174 {
 175     int i, offset;
 176     bool is_periodic = false;
 177
 178     if (value_notzero_p(e->d))
 179         return false;
 180
 181     assert(e->x.p->type != partition);
 182     if (e->x.p->type == periodic) {
 183         Value size;
 184         value_init(size);
 185         value_set_si(size, e->x.p->size);
 186         value_lcm(periods->p[e->x.p->pos-1], periods->p[e->x.p->pos-1], size);
 187         value_clear(size);
 188         is_periodic = true;
 189     }
 190     offset = type_offset(e->x.p);
 191     for (i = e->x.p->size-1; i >= offset; --i)
 192         is_periodic = evalue_is_periodic(&e->x.p->arr[i], periods) || is_periodic;
 193     return is_periodic;
 194 }
 195
 196 static ex evalue2lst(const evalue *e, const exvector& vars,
 197                      exvector& extravar, vector<typed_evalue>& expr,
 198                      Vector *periods)
 199 {
 200     Vector *coset = Vector_Alloc(periods->Size);
 201     lst list;
 202     for (;;) {
 203         int i;
 204         list.append(evalue2ex_r(e, vars, extravar, expr, coset));
 205         for (i = coset->Size-1; i >= 0; --i) {
 206             value_increment(coset->p[i], coset->p[i]);
 207             if (value_lt(coset->p[i], periods->p[i]))
 208                 break;
 209             value_set_si(coset->p[i], 0);
 210         }
 211         if (i < 0)
 212             break;
 213     }
 214     Vector_Free(coset);
 215     return list;
 216 }
 217
 218 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar,
 219              Matrix **C, Vector **p)
 220 {
 221     vector<typed_evalue> expr;
 222     Vector *periods = Vector_Alloc(vars.size());
 223     assert(p);
 224     assert(C);
 225     for (int i = 0; i < periods->Size; ++i)
 226         value_set_si(periods->p[i], 1);
 227     if (evalue_is_periodic(e, periods)) {
 228         *p = periods;
 229         *C = NULL;
 230         lst list;
 231         return list;
 232     } else {
 233         Vector_Free(periods);
 234         *p = NULL;
 235         ex poly = evalue2ex_r(e, vars, floorvar, expr, NULL);
 236         Matrix *M = setup_constraints(expr, vars.size());
 237         *C = M;
 238         return poly;
 239     }
 240 }
 241
 242 /* if the evalue is a relation, we use the relation to cut off the
 243  * the edges of the domain
 244  */
 245 static Polyhedron *relation_domain(Polyhedron *D, evalue *fr, unsigned MaxRays)
 246 {
 247     assert(value_zero_p(fr->d));
 248     assert(fr->x.p->type == fractional);
 249     assert(fr->x.p->size == 3);
 250     Matrix *T = Matrix_Alloc(2, D->Dimension+1);
 251     value_set_si(T->p[1][D->Dimension], 1);
 252
 253     /* convert argument of fractional to polylib */
 254     /* the argument is assumed to be linear */
 255     evalue *p = &fr->x.p->arr[0];
 256     evalue_denom(p, &T->p[1][D->Dimension]);
 257     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
 258         assert(p->x.p->type == polynomial);
 259         assert(p->x.p->size == 2);
 260         assert(value_notzero_p(p->x.p->arr[1].d));
 261         int pos = p->x.p->pos - 1;
 262         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
 263         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][D->Dimension]);
 264         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
 265     }
 266     int pos = D->Dimension;
 267     value_assign(T->p[0][pos], p->x.n);
 268     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][D->Dimension]);
 269     value_division(T->p[0][pos], T->p[0][pos], p->d);
 270
 271     Polyhedron *E = NULL;
 272     for (Polyhedron *P = D; P; P = P->next) {
 273         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
 274         I = DomainConstraintSimplify(I, MaxRays);
 275         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
 276         Polyhedron_Free(I);
 277         Polyhedron *next = P->next;
 278         P->next = NULL;
 279         Polyhedron *S = DomainIntersection(P, R, MaxRays);
 280         Polyhedron_Free(R);
 281         P->next = next;
 282         if (emptyQ2(S))
 283             Polyhedron_Free(S);
 284         else
 285             E = DomainConcat(S, E);
 286     }
 287     Matrix_Free(T);
 288
 289     return E;
 290 }
 291
 292 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 293                                       Polyhedron *ctx, const exvector& params)
 294 {
 295     piecewise_lst *pl;
 296     barvinok_options *options = barvinok_options_new_with_defaults();
 297     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
 298     barvinok_options_free(options);
 299     return pl;
 300 }
 301
 302 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 303                             Polyhedron *D, const ex& poly,
 304                             Polyhedron *ctx,
 305                             const exvector& params, const exvector& floorvar,
 306                             barvinok_options *options);
 307
 308 /* Recursively apply Bernstein expansion on P, optimizing over dims[i]
 309  * variables in each level.  The context ctx is assumed to have been adapted
 310  * to the first level in the recursion.
 311  */
 312 static piecewise_lst *bernstein_coefficients_recursive(piecewise_lst *pl_all,
 313                             Polyhedron *P, const vector<int>& dims, const ex& poly,
 314                             Polyhedron *ctx,
 315                             const exvector& params, const exvector& vars,
 316                             barvinok_options *options)
 317 {
 318     assert(dims.size() > 0);
 319     assert(ctx->Dimension == P->Dimension - dims[0]);
 320     piecewise_lst *pl;
 321     unsigned done = 0;
 322     for (int j = 0; j < dims.size(); ++j) {
 323         exvector pl_vars;
 324         pl_vars.insert(pl_vars.end(), vars.begin()+done, vars.begin()+done+dims[j]);
 325         exvector pl_params;
 326         pl_params.insert(pl_params.end(), vars.begin()+done+dims[j], vars.end());
 327         pl_params.insert(pl_params.end(), params.begin(), params.end());
 328
 329         if (!j)
 330             pl = bernstein_coefficients(NULL, P, poly, ctx,
 331                                         pl_params, pl_vars, options);
 332         else {
 333             piecewise_lst *new_pl = NULL;
 334             Polyhedron *U = Universe_Polyhedron(pl_params.size());
 335
 336             for (int i = 0; i < pl->list.size(); ++i) {
 337                 Polyhedron *D = pl->list[i].first;
 338                 lst polys = pl->list[i].second;
 339                 new_pl = bernstein_coefficients(new_pl, D, polys, U, pl_params,
 340                                                 pl_vars, options);
 341             }
 342
 343             Polyhedron_Free(U);
 344
 345             delete pl;
 346             pl = new_pl;
 347         }
 348
 349         done += dims[j];
 350
 351         if (!pl)
 352             return pl_all;
 353     }
 354
 355     if (!pl_all)
 356         pl_all = pl;
 357     else {
 358         pl_all->combine(*pl);
 359         delete pl;
 360     }
 361
 362     return pl_all;
 363 }
 364
 365 static piecewise_lst *bernstein_coefficients_full_recurse(piecewise_lst *pl_all,
 366                             Polyhedron *P, const ex& poly,
 367                             Polyhedron *ctx,
 368                             const exvector& params, const exvector& vars,
 369                             barvinok_options *options)
 370 {
 371     Polyhedron *CR = align_context(ctx, P->Dimension-1, options->MaxRays);
 372     vector<int> dims(vars.size());
 373     for (int i = 0; i < dims.size(); ++i)
 374         dims[i] = 1;
 375     pl_all = bernstein_coefficients_recursive(pl_all, P, dims, poly, CR,
 376                                               params, vars, options);
 377     Polyhedron_Free(CR);
 378
 379     return pl_all;
 380 }
 381
 382 static piecewise_lst *bernstein_coefficients_product(piecewise_lst *pl_all,
 383                             Polyhedron *F, Matrix *T, const ex& poly,
 384                             Polyhedron *ctx,
 385                             const exvector& params, const exvector& vars,
 386                             barvinok_options *options)
 387 {
 388     if (emptyQ2(ctx))
 389         return pl_all;
 390     for (Polyhedron *G = F; G; G = G->next)
 391         if (emptyQ2(G))
 392             return pl_all;
 393
 394     unsigned nparam = params.size();
 395     unsigned nvar = vars.size();
 396     unsigned constraints;
 397     unsigned factors;
 398     Polyhedron *C = NULL;
 399
 400     /* More context constraints */
 401     if (F->Dimension == ctx->Dimension) {
 402         C = F;
 403         F = F->next;
 404     }
 405     assert(F);
 406     assert(F->next);
 407
 408     Matrix *M;
 409     Polyhedron *P;
 410     Polyhedron *PC;
 411     M = Matrix_Alloc(F->NbConstraints, 1+nvar+nparam+1);
 412     for (int i = 0; i < F->NbConstraints; ++i) {
 413         Vector_Copy(F->Constraint[i], M->p[i], 1+F->Dimension-nparam);
 414         Vector_Copy(F->Constraint[i]+1+F->Dimension-nparam,
 415                     M->p[i]+1+nvar, nparam+1);
 416     }
 417     P = Constraints2Polyhedron(M, options->MaxRays);
 418     Matrix_Free(M);
 419
 420     factors = 1;
 421     constraints = C ? C->NbConstraints : 0;
 422     constraints += ctx->NbConstraints;
 423     for (Polyhedron *G = F->next; G; G = G->next) {
 424         constraints += G->NbConstraints;
 425         ++factors;
 426     }
 427
 428     unsigned total_var = nvar-(F->Dimension-nparam);
 429     unsigned skip = 0;
 430     unsigned c = 0;
 431     M = Matrix_Alloc(constraints, 1+total_var+nparam+1);
 432     for (Polyhedron *G = F->next; G; G = G->next) {
 433         unsigned this_var = G->Dimension - nparam;
 434         for (int i = 0; i < G->NbConstraints; ++i) {
 435             value_assign(M->p[c+i][0], G->Constraint[i][0]);
 436             Vector_Copy(G->Constraint[i]+1, M->p[c+i]+1+skip, this_var);
 437             Vector_Copy(G->Constraint[i]+1+this_var, M->p[c+i]+1+total_var,
 438                         nparam+1);
 439         }
 440         c += G->NbConstraints;
 441         skip += this_var;
 442     }
 443     assert(skip == total_var);
 444     if (C) {
 445         for (int i = 0; i < C->NbConstraints; ++i) {
 446             value_assign(M->p[c+i][0], C->Constraint[i][0]);
 447             Vector_Copy(C->Constraint[i]+1, M->p[c+i]+1+total_var,
 448                         nparam+1);
 449         }
 450         c += C->NbConstraints;
 451     }
 452     for (int i = 0; i < ctx->NbConstraints; ++i) {
 453         value_assign(M->p[c+i][0], ctx->Constraint[i][0]);
 454         Vector_Copy(ctx->Constraint[i]+1, M->p[c+i]+1+total_var, nparam+1);
 455     }
 456     PC = Constraints2Polyhedron(M, options->MaxRays);
 457     Matrix_Free(M);
 458
 459     exvector newvars = constructVariableVector(nvar, "t");
 460     matrix subs(1, nvar);
 461     for (int i = 0; i < nvar; ++i)
 462         for (int j = 0; j < nvar; ++j)
 463             subs(0,i) += value2numeric(T->p[i][j]) * newvars[j];
 464
 465     ex newpoly = replaceVariablesInPolynomial(poly, vars, subs);
 466
 467     vector<int> dims(factors);
 468     for (int i = 0; F; ++i, F = F->next)
 469         dims[i] = F->Dimension-nparam;
 470
 471     pl_all = bernstein_coefficients_recursive(pl_all, P, dims, newpoly, PC,
 472                                               params, newvars, options);
 473
 474     Polyhedron_Free(P);
 475     Polyhedron_Free(PC);
 476
 477     return pl_all;
 478 }
 479
 480 static piecewise_lst *bernstein_coefficients_polyhedron(piecewise_lst *pl_all,
 481                             Polyhedron *P, const ex& poly,
 482                             Polyhedron *ctx,
 483                             const exvector& params, const exvector& floorvar,
 484                             barvinok_options *options)
 485 {
 486     if (Polyhedron_is_unbounded(P, ctx->Dimension, options->MaxRays)) {
 487         fprintf(stderr, "warning: unbounded domain skipped\n");
 488         Polyhedron_Print(stderr, P_VALUE_FMT, P);
 489         return pl_all;
 490     }
 491
 492     if (options->bernstein_recurse & BV_BERNSTEIN_FACTORS) {
 493         Matrix *T = NULL;
 494         Polyhedron *F = Polyhedron_Factor(P, ctx->Dimension, &T, options->MaxRays);
 495         if (F) {
 496             pl_all = bernstein_coefficients_product(pl_all, F, T, poly, ctx, params,
 497                                                     floorvar, options);
 498             Domain_Free(F);
 499             Matrix_Free(T);
 500             return pl_all;
 501         }
 502     }
 503     if (floorvar.size() > 1 &&
 504                 options->bernstein_recurse & BV_BERNSTEIN_INTERVALS)
 505         return bernstein_coefficients_full_recurse(pl_all, P, poly, ctx,
 506                                                    params, floorvar, options);
 507
 508     unsigned PP_MaxRays = options->MaxRays;
 509     if (PP_MaxRays & POL_NO_DUAL)
 510         PP_MaxRays = 0;
 511
 512     Param_Polyhedron *PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 513     if (!PP)
 514         return pl_all;
 515     piecewise_lst *pl = new piecewise_lst(params, options->bernstein_optimize);
 516
 517     int nd = -1;
 518     Polyhedron *TC = true_context(P, ctx, options->MaxRays);
 519     FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, PD, rVD)
 520         matrix VM = domainVertices(PP, PD, params);
 521         lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 522         pl->add_guarded_lst(rVD, coeffs);
 523     END_FORALL_REDUCED_DOMAIN
 524     Polyhedron_Free(TC);
 525
 526     Param_Polyhedron_Free(PP);
 527     if (!pl_all)
 528         pl_all = pl;
 529     else {
 530         pl_all->combine(*pl);
 531         delete pl;
 532     }
 533
 534     return pl_all;
 535 }
 536
 537 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 538                             Polyhedron *D, const ex& poly,
 539                             Polyhedron *ctx,
 540                             const exvector& params, const exvector& floorvar,
 541                             barvinok_options *options)
 542 {
 543     if (!D->next && emptyQ2(D))
 544         return pl_all;
 545
 546     for (Polyhedron *P = D; P; P = P->next) {
 547         /* This shouldn't happen */
 548         if (emptyQ2(P))
 549             continue;
 550         Polyhedron *next = P->next;
 551         P->next = NULL;
 552         pl_all = bernstein_coefficients_polyhedron(pl_all, P, poly, ctx,
 553                                                    params, floorvar, options);
 554         P->next = next;
 555     }
 556     return pl_all;
 557 }
 558
 559 /* Compute the coefficients of the polynomial corresponding to each coset
 560  * on its own domain.  This allows us to cut the domain on multiples of
 561  * the period.
 562  * To perform the cutting for a coset "i mod n = c" we map the domain
 563  * to the quotient space trough "i = i' n + c", simplify the constraints
 564  * (implicitly) and then map back to the original space.
 565  */
 566 static piecewise_lst *bernstein_coefficients_periodic(piecewise_lst *pl_all,
 567                             Polyhedron *D, const evalue *e,
 568                             Polyhedron *ctx, const exvector& vars,
 569                             const exvector& params, Vector *periods,
 570                             barvinok_options *options)
 571 {
 572     assert(D->Dimension == periods->Size);
 573     Matrix *T = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 574     Matrix *T2 = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 575     Vector *coset = Vector_Alloc(periods->Size);
 576     exvector extravar;
 577     vector<typed_evalue> expr;
 578     exvector allvars = vars;
 579     allvars.insert(allvars.end(), params.begin(), params.end());
 580
 581     value_set_si(T2->p[D->Dimension][D->Dimension], 1);
 582     for (int i = 0; i < D->Dimension; ++i) {
 583         value_assign(T->p[i][i], periods->p[i]);
 584         value_lcm(T2->p[D->Dimension][D->Dimension],
 585                   T2->p[D->Dimension][D->Dimension], periods->p[i]);
 586     }
 587     value_set_si(T->p[D->Dimension][D->Dimension], 1);
 588     for (int i = 0; i < D->Dimension; ++i)
 589         value_division(T2->p[i][i], T2->p[D->Dimension][D->Dimension],
 590                        periods->p[i]);
 591     for (;;) {
 592         int i;
 593         ex poly = evalue2ex_r(e, allvars, extravar, expr, coset);
 594         assert(extravar.size() == 0);
 595         assert(expr.size() == 0);
 596         Polyhedron *E = DomainPreimage(D, T, options->MaxRays);
 597         Polyhedron *F = DomainPreimage(E, T2, options->MaxRays);
 598         Polyhedron_Free(E);
 599         if (!emptyQ(F))
 600             pl_all = bernstein_coefficients(pl_all, F, poly, ctx, params,
 601                                             vars, options);
 602         Polyhedron_Free(F);
 603         for (i = D->Dimension-1; i >= 0; --i) {
 604             value_increment(coset->p[i], coset->p[i]);
 605             value_increment(T->p[i][D->Dimension], T->p[i][D->Dimension]);
 606             value_subtract(T2->p[i][D->Dimension], T2->p[i][D->Dimension],
 607                            T2->p[i][i]);
 608             if (value_lt(coset->p[i], periods->p[i]))
 609                 break;
 610             value_set_si(coset->p[i], 0);
 611             value_set_si(T->p[i][D->Dimension], 0);
 612             value_set_si(T2->p[i][D->Dimension], 0);
 613         }
 614         if (i < 0)
 615             break;
 616     }
 617     Vector_Free(coset);
 618     Matrix_Free(T);
 619     Matrix_Free(T2);
 620     return pl_all;
 621 }
 622
 623 piecewise_lst *bernstein_coefficients_relation(piecewise_lst *pl_all,
 624                     Polyhedron *D, evalue *EP, Polyhedron *ctx,
 625                     const exvector& allvars, const exvector& vars,
 626                     const exvector& params, barvinok_options *options)
 627 {
 628     if (value_zero_p(EP->d) && EP->x.p->type == relation) {
 629         Polyhedron *E = relation_domain(D, &EP->x.p->arr[0], options->MaxRays);
 630         if (E) {
 631             pl_all = bernstein_coefficients_relation(pl_all, E, &EP->x.p->arr[1],
 632                                                      ctx, allvars, vars, params,
 633                                                      options);
 634             Domain_Free(E);
 635         }
 636         /* In principle, we could cut off the edges of this domain too */
 637         if (EP->x.p->size > 2)
 638             pl_all = bernstein_coefficients_relation(pl_all, D, &EP->x.p->arr[2],
 639                                                      ctx, allvars, vars, params,
 640                                                      options);
 641         return pl_all;
 642     }
 643
 644     Matrix *M;
 645     exvector floorvar;
 646     Vector *periods;
 647     ex poly = evalue2ex(EP, allvars, floorvar, &M, &periods);
 648     floorvar.insert(floorvar.end(), vars.begin(), vars.end());
 649     Polyhedron *E = D;
 650     if (M) {
 651         Polyhedron *AE = align_context(D, M->NbColumns-2, options->MaxRays);
 652         E = DomainAddConstraints(AE, M, options->MaxRays);
 653         Matrix_Free(M);
 654         Domain_Free(AE);
 655     }
 656     if (is_exactly_a<fail>(poly)) {
 657         delete pl_all;
 658         return NULL;
 659     }
 660     if (periods)
 661         pl_all = bernstein_coefficients_periodic(pl_all, E, EP, ctx, vars,
 662                                                  params, periods, options);
 663     else
 664         pl_all = bernstein_coefficients(pl_all, E, poly, ctx, params,
 665                                         floorvar, options);
 666     if (periods)
 667         Vector_Free(periods);
 668     if (D != E)
 669         Domain_Free(E);
 670
 671     return pl_all;
 672 }
 673
 674 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 675                                       Polyhedron *ctx, const exvector& params,
 676                                       barvinok_options *options)
 677 {
 678     unsigned nparam = ctx->Dimension;
 679     if (EVALUE_IS_ZERO(*e))
 680         return pl_all;
 681     assert(value_zero_p(e->d));
 682     assert(e->x.p->type == partition);
 683     assert(e->x.p->size >= 2);
 684     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 685
 686     exvector vars = constructVariableVector(nvars, "v");
 687     exvector allvars = vars;
 688     allvars.insert(allvars.end(), params.begin(), params.end());
 689
 690     for (int i = 0; i < e->x.p->size/2; ++i) {
 691         pl_all = bernstein_coefficients_relation(pl_all,
 692                         EVALUE_DOMAIN(e->x.p->arr[2*i]), &e->x.p->arr[2*i+1],
 693                         ctx, allvars, vars, params, options);
 694     }
 695     return pl_all;
 696 }
 697
 698 static __isl_give isl_qpolynomial *qp_from_ex(__isl_take isl_dim *dim,
 699         const GiNaC::ex ex, const GiNaC::exvector &params, int i)
 700 {
 701         isl_qpolynomial *qp;
 702         isl_qpolynomial *base;
 703         int deg;
 704         int j;
 705
 706         if (is_a<fail>(ex))
 707                 return isl_qpolynomial_nan(dim);
 708
 709         if (is_a<numeric>(ex)) {
 710                 numeric r = ex_to<numeric>(ex);
 711                 isl_int n;
 712                 isl_int d;
 713                 isl_int_init(n);
 714                 isl_int_init(d);
 715                 numeric2value(r.numer(), n);
 716                 numeric2value(r.denom(), d);
 717                 qp = isl_qpolynomial_rat_cst(dim, n, d);
 718                 isl_int_clear(n);
 719                 isl_int_clear(d);
 720                 return qp;
 721         }
 722
 723         deg = ex.degree(params[i]);
 724         if (deg == 0)
 725                 return qp_from_ex(dim, ex, params, i + 1);
 726
 727         base = isl_qpolynomial_var(isl_dim_copy(dim), isl_dim_param, i);
 728         qp = qp_from_ex(isl_dim_copy(dim), ex.coeff(params[i], deg),
 729                         params, i + 1);
 730
 731         for (j = deg - 1; j >= 0; --j) {
 732                 qp = isl_qpolynomial_mul(qp, isl_qpolynomial_copy(base));
 733                 qp = isl_qpolynomial_add(qp,
 734                             qp_from_ex(isl_dim_copy(dim), ex.coeff(params[i], j),
 735                                         params, i + 1));
 736         }
 737
 738         isl_qpolynomial_free(base);
 739         isl_dim_free(dim);
 740
 741         return qp;
 742 }
 743
 744 __isl_give isl_qpolynomial *isl_qpolynomial_from_ginac(__isl_take isl_dim *dim,
 745         const GiNaC::ex &ex, const GiNaC::exvector &params)
 746 {
 747         GiNaC::ex exp;
 748
 749         if (!dim)
 750                 return NULL;
 751
 752         exp = ex.expand();
 753         return qp_from_ex(dim, exp, params, 0);
 754 error:
 755         isl_dim_free(dim);
 756         return NULL;
 757 }
 758
 759 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_from_ginac(
 760         __isl_take isl_dim *dim, enum isl_fold type, const GiNaC::lst &lst,
 761         const GiNaC::exvector &params)
 762 {
 763         isl_qpolynomial_fold *fold;
 764         lst::const_iterator j;
 765
 766         fold = isl_qpolynomial_fold_empty(type, isl_dim_copy(dim));
 767         for (j = lst.begin(); j != lst.end(); ++j) {
 768                 isl_qpolynomial *qp;
 769                 isl_qpolynomial_fold *fold_i;
 770
 771                 qp = isl_qpolynomial_from_ginac(isl_dim_copy(dim), *j, params);
 772                 fold_i = isl_qpolynomial_fold_alloc(type, qp);
 773                 fold = isl_qpolynomial_fold_fold(fold, fold_i);
 774         }
 775         isl_dim_free(dim);
 776         return fold;
 777 }
 778
 779 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_ginac(
 780         __isl_take isl_dim *dim, bernstein::piecewise_lst *pl,
 781         const GiNaC::exvector &params)
 782 {
 783         enum isl_fold type;
 784         isl_pw_qpolynomial_fold *pwf;
 785
 786         pwf = isl_pw_qpolynomial_fold_zero(isl_dim_copy(dim));
 787
 788         type = pl->sign > 0 ? isl_fold_max
 789                             : pl->sign < 0 ? isl_fold_min : isl_fold_list;
 790
 791         for (int i = 0; i < pl->list.size(); ++i) {
 792                 isl_pw_qpolynomial_fold *pwf_i;
 793                 isl_qpolynomial_fold *fold;
 794                 isl_set *set;
 795
 796                 set = isl_set_new_from_polylib(pl->list[i].first,
 797                                 isl_dim_copy(dim));
 798                 fold = isl_qpolynomial_fold_from_ginac(isl_dim_copy(dim),
 799                                 type, pl->list[i].second, params);
 800                 pwf_i = isl_pw_qpolynomial_fold_alloc(set, fold);
 801                 pwf = isl_pw_qpolynomial_fold_add_disjoint(pwf, pwf_i);
 802         }
 803
 804         isl_dim_free(dim);
 805
 806         return pwf;
 807 }
 808
 809 }
 810
 811 struct isl_bound {
 812         piecewise_lst *pl;
 813         Polyhedron *U;
 814         exvector vars;
 815         exvector params;
 816         exvector allvars;
 817 };
 818
 819 static int guarded_qp_bernstein_coefficients(__isl_take isl_set *set,
 820         __isl_take isl_qpolynomial *qp, void *user)
 821 {
 822         struct isl_bound *bound = (struct isl_bound *)user;
 823         unsigned nvar;
 824         Polyhedron *D;
 825         struct barvinok_options *options;
 826         evalue *e;
 827
 828         options = barvinok_options_new_with_defaults();
 829
 830         if (!set || !qp)
 831                 goto error;
 832
 833         nvar = isl_set_dim(set, isl_dim_set);
 834
 835         e = isl_qpolynomial_to_evalue(qp);
 836         if (!e)
 837                 goto error;
 838
 839         set = isl_set_make_disjoint(set);
 840         D = isl_set_to_polylib(set);
 841
 842         bound->vars = constructVariableVector(nvar, "v");
 843         bound->allvars = bound->vars;
 844         bound->allvars.insert(bound->allvars.end(),
 845                                 bound->params.begin(), bound->params.end());
 846
 847         bound->pl = bernstein_coefficients_relation(bound->pl, D, e, bound->U,
 848                         bound->allvars, bound->vars, bound->params, options);
 849
 850         Domain_Free(D);
 851         evalue_free(e);
 852         isl_set_free(set);
 853         isl_qpolynomial_free(qp);
 854         barvinok_options_free(options);
 855
 856         return 0;
 857 error:
 858         isl_set_free(set);
 859         isl_qpolynomial_free(qp);
 860         barvinok_options_free(options);
 861         return -1;
 862 }
 863
 864 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_bound_bernstein(
 865         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type)
 866 {
 867         isl_dim *dim = NULL;
 868         unsigned nvar;
 869         unsigned nparam;
 870         struct isl_bound bound;
 871         struct isl_pw_qpolynomial_fold *pwf;
 872
 873         if (!pwqp)
 874                 return NULL;
 875
 876         dim = isl_pw_qpolynomial_get_dim(pwqp);
 877         nvar = isl_dim_size(dim, isl_dim_set);
 878
 879         if (isl_pw_qpolynomial_is_zero(pwqp)) {
 880                 isl_pw_qpolynomial_free(pwqp);
 881                 dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 882                 return isl_pw_qpolynomial_fold_zero(dim);
 883         }
 884         if (nvar == 0) {
 885                 isl_dim_free(dim);
 886                 return isl_pw_qpolynomial_fold_from_pw_qpolynomial(type, pwqp);
 887         }
 888
 889         nparam = isl_dim_size(dim, isl_dim_param);
 890         bound.pl = NULL;
 891         bound.U = Universe_Polyhedron(nparam);
 892         bound.params = constructVariableVector(nparam, "p");
 893
 894         if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp,
 895                                     guarded_qp_bernstein_coefficients, &bound))
 896                 goto error;
 897
 898         if (type == isl_fold_max)
 899                 bound.pl->maximize();
 900         else
 901                 bound.pl->minimize();
 902
 903         dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 904
 905         if (type == isl_fold_max)
 906                 bound.pl->sign = 1;
 907         else
 908                 bound.pl->sign = -1;
 909         pwf = isl_pw_qpolynomial_fold_from_ginac(dim, bound.pl, bound.params);
 910
 911         Polyhedron_Free(bound.U);
 912         delete bound.pl;
 913
 914         isl_pw_qpolynomial_free(pwqp);
 915
 916         return pwf;
 917 error:
 918         isl_pw_qpolynomial_free(pwqp);
 919         isl_dim_free(dim);
 920         return NULL;
 921 }