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