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