bernstein.cc

   1 #include <bernstein/bernstein.h>
   2 #include <bernstein/piecewise_lst.h>
   3 #include <barvinok/barvinok.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/bernstein.h>
   6 #include <barvinok/options.h>
   7
   8 using namespace GiNaC;
   9 using namespace bernstein;
  10
  11 namespace barvinok {
  12
  13 ex evalue2ex(evalue *e, const exvector& vars)
  14 {
  15     if (value_notzero_p(e->d))
  16         return value2numeric(e->x.n)/value2numeric(e->d);
  17     if (e->x.p->type != polynomial)
  18         return fail();
  19     ex poly = 0;
  20     for (int i = e->x.p->size-1; i >= 0; --i) {
  21         poly *= vars[e->x.p->pos-1];
  22         ex t = evalue2ex(&e->x.p->arr[i], vars);
  23         if (is_exactly_a<fail>(t))
  24             return t;
  25         poly += t;
  26     }
  27     return poly;
  28 }
  29
  30 /* if the evalue is a relation, we use the relation to cut off the
  31  * the edges of the domain
  32  */
  33 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
  34                                 unsigned MaxRays)
  35 {
  36     *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
  37     *poly = e = &e->x.p->arr[2*i+1];
  38     if (value_notzero_p(e->d))
  39         return;
  40     if (e->x.p->type != relation)
  41         return;
  42     if (e->x.p->size > 2)
  43         return;
  44     evalue *fr = &e->x.p->arr[0];
  45     assert(value_zero_p(fr->d));
  46     assert(fr->x.p->type == fractional);
  47     assert(fr->x.p->size == 3);
  48     Matrix *T = Matrix_Alloc(2, (*D)->Dimension+1);
  49     value_set_si(T->p[1][(*D)->Dimension], 1);
  50
  51     /* convert argument of fractional to polylib */
  52     /* the argument is assumed to be linear */
  53     evalue *p = &fr->x.p->arr[0];
  54     evalue_denom(p, &T->p[1][(*D)->Dimension]);
  55     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
  56         assert(p->x.p->type == polynomial);
  57         assert(p->x.p->size == 2);
  58         assert(value_notzero_p(p->x.p->arr[1].d));
  59         int pos = p->x.p->pos - 1;
  60         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
  61         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
  62         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
  63     }
  64     int pos = (*D)->Dimension;
  65     value_assign(T->p[0][pos], p->x.n);
  66     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
  67     value_division(T->p[0][pos], T->p[0][pos], p->d);
  68
  69     Polyhedron *E = NULL;
  70     for (Polyhedron *P = *D; P; P = P->next) {
  71         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
  72         I = DomainConstraintSimplify(I, MaxRays);
  73         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
  74         Polyhedron_Free(I);
  75         Polyhedron *next = P->next;
  76         P->next = NULL;
  77         Polyhedron *S = DomainIntersection(P, R, MaxRays);
  78         Polyhedron_Free(R);
  79         P->next = next;
  80         if (emptyQ2(S))
  81             Polyhedron_Free(S);
  82         else
  83             E = DomainConcat(S, E);
  84     }
  85     Matrix_Free(T);
  86
  87     *D = E;
  88     *poly = &e->x.p->arr[1];
  89 }
  90
  91 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
  92                                       Polyhedron *ctx, const exvector& params)
  93 {
  94     piecewise_lst *pl;
  95     barvinok_options *options = barvinok_options_new_with_defaults();
  96     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
  97     barvinok_options_free(options);
  98     return pl;
  99 }
 100
 101 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 102                                       Polyhedron *ctx, const exvector& params,
 103                                       barvinok_options *options)
 104 {
 105     unsigned nparam = ctx->Dimension;
 106     if (EVALUE_IS_ZERO(*e))
 107         return pl_all;
 108     assert(value_zero_p(e->d));
 109     assert(e->x.p->type == partition);
 110     assert(e->x.p->size >= 2);
 111     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 112     unsigned PP_MaxRays = options->MaxRays;
 113     if (PP_MaxRays & POL_NO_DUAL)
 114         PP_MaxRays = 0;
 115
 116     exvector vars = constructVariableVector(nvars, "v");
 117     exvector allvars = vars;
 118     allvars.insert(allvars.end(), params.begin(), params.end());
 119
 120     for (int i = 0; i < e->x.p->size/2; ++i) {
 121         Param_Polyhedron *PP;
 122         Polyhedron *E;
 123         evalue *EP;
 124         evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
 125         ex poly = evalue2ex(EP, allvars);
 126         if (is_exactly_a<fail>(poly)) {
 127             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 128                 Domain_Free(E);
 129             delete pl_all;
 130             return NULL;
 131         }
 132         for (Polyhedron *P = E; P; P = P->next) {
 133             Polyhedron *next = P->next;
 134             piecewise_lst *pl = new piecewise_lst(params);
 135             Polyhedron *P1 = P;
 136             P->next = NULL;
 137             PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 138             for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 139                 matrix VM = domainVertices(PP, Q, params);
 140                 lst coeffs = bernsteinExpansion(VM, poly, vars, params);
 141                 pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
 142             }
 143             Param_Polyhedron_Free(PP);
 144             if (!pl_all)
 145                 pl_all = pl;
 146             else {
 147                 pl_all->combine(*pl);
 148                 delete pl;
 149             }
 150             P->next = next;
 151         }
 152         if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 153             Domain_Free(E);
 154     }
 155     return pl_all;
 156 }
 157
 158 }