bernstein.cc

   1 #include <vector>
   2 #include <bernstein/bernstein.h>
   3 #include <bernstein/piecewise_lst.h>
   4 #include <barvinok/barvinok.h>
   5 #include <barvinok/util.h>
   6 #include <barvinok/bernstein.h>
   7 #include <barvinok/options.h>
   8
   9 using namespace GiNaC;
  10 using namespace bernstein;
  11
  12 using std::pair;
  13 using std::vector;
  14 using std::cerr;
  15 using std::endl;
  16
  17 namespace barvinok {
  18
  19 ex evalue2ex(evalue *e, const exvector& vars)
  20 {
  21     if (value_notzero_p(e->d))
  22         return value2numeric(e->x.n)/value2numeric(e->d);
  23     if (e->x.p->type != polynomial)
  24         return fail();
  25     ex poly = 0;
  26     for (int i = e->x.p->size-1; i >= 0; --i) {
  27         poly *= vars[e->x.p->pos-1];
  28         ex t = evalue2ex(&e->x.p->arr[i], vars);
  29         if (is_exactly_a<fail>(t))
  30             return t;
  31         poly += t;
  32     }
  33     return poly;
  34 }
  35
  36 static int type_offset(enode *p)
  37 {
  38    return p->type == fractional ? 1 :
  39           p->type == flooring ? 1 : 0;
  40 }
  41
  42 typedef pair<bool, const evalue *> typed_evalue;
  43
  44 static ex evalue2ex_add_var(evalue *e, exvector& extravar,
  45                             vector<typed_evalue>& expr, bool is_fract)
  46 {
  47     ex base_var = 0;
  48
  49     for (int i = 0; i < expr.size(); ++i) {
  50         if (is_fract == expr[i].first && eequal(e, expr[i].second)) {
  51             base_var = extravar[i];
  52             break;
  53         }
  54     }
  55     if (base_var != 0)
  56         return base_var;
  57
  58     char name[20];
  59     snprintf(name, sizeof(name), "f%c%d", is_fract ? 'r' : 'l', expr.size());
  60     extravar.push_back(base_var = symbol(name));
  61     expr.push_back(typed_evalue(is_fract, e));
  62
  63     return base_var;
  64 }
  65
  66 static ex evalue2ex_r(const evalue *e, const exvector& vars,
  67                       exvector& extravar, vector<typed_evalue>& expr)
  68 {
  69     if (value_notzero_p(e->d))
  70         return value2numeric(e->x.n)/value2numeric(e->d);
  71     ex base_var = 0;
  72     ex poly = 0;
  73
  74     switch (e->x.p->type) {
  75     case polynomial:
  76         base_var = vars[e->x.p->pos-1];
  77         break;
  78     case flooring:
  79         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, false);
  80         break;
  81     case fractional:
  82         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, true);
  83         break;
  84     default:
  85         return fail();
  86     }
  87
  88     int offset = type_offset(e->x.p);
  89     for (int i = e->x.p->size-1; i >= offset; --i) {
  90         poly *= base_var;
  91         ex t = evalue2ex_r(&e->x.p->arr[i], vars, extravar, expr);
  92         if (is_exactly_a<fail>(t))
  93             return t;
  94         poly += t;
  95     }
  96     return poly;
  97 }
  98
  99 /* For each t = floor(e/d), set up two constraints
 100  *
 101  *       e - d t       >= 0
 102  *      -e + d t + d-1 >= 0
 103  *
 104  * e is assumed to be an affine expression.
 105  *
 106  * For each t = fract(e/d), set up two constraints
 107  *
 108  *          -d t + d-1 >= 0
 109  *             t       >= 0
 110  */
 111 static Matrix *setup_constraints(const vector<typed_evalue> expr, int nvar)
 112 {
 113     int extra = expr.size();
 114     if (!extra)
 115         return NULL;
 116     Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
 117     for (int i = 0; i < extra; ++i) {
 118         Value *d = &M->p[2*i][1+i];
 119         value_set_si(*d, 1);
 120         evalue_denom(expr[i].second, d);
 121         if (expr[i].first) {
 122             value_set_si(M->p[2*i][0], 1);
 123             value_decrement(M->p[2*i][1+extra+nvar], *d);
 124             value_oppose(*d, *d);
 125             value_set_si(M->p[2*i+1][0], 1);
 126             value_set_si(M->p[2*i+1][1+i], 1);
 127         } else {
 128             const evalue *e;
 129             for (e = expr[i].second; value_zero_p(e->d); e = &e->x.p->arr[0]) {
 130                 assert(e->x.p->type == polynomial);
 131                 assert(e->x.p->size == 2);
 132                 evalue *c = &e->x.p->arr[1];
 133                 value_multiply(M->p[2*i][1+extra+e->x.p->pos-1], *d, c->x.n);
 134                 value_division(M->p[2*i][1+extra+e->x.p->pos-1],
 135                                M->p[2*i][1+extra+e->x.p->pos-1], c->d);
 136             }
 137             value_multiply(M->p[2*i][1+extra+nvar], *d, e->x.n);
 138             value_division(M->p[2*i][1+extra+nvar], M->p[2*i][1+extra+nvar], e->d);
 139             value_oppose(*d, *d);
 140             value_set_si(M->p[2*i][0], -1);
 141             Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
 142             value_set_si(M->p[2*i][0], 1);
 143             value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
 144             value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
 145         }
 146     }
 147     return M;
 148 }
 149
 150 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar, Matrix **C)
 151 {
 152     vector<typed_evalue> expr;
 153     ex poly = evalue2ex_r(e, vars, floorvar, expr);
 154     assert(C);
 155     Matrix *M = setup_constraints(expr, vars.size());
 156     *C = M;
 157     return poly;
 158 }
 159
 160 /* if the evalue is a relation, we use the relation to cut off the
 161  * the edges of the domain
 162  */
 163 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
 164                                 unsigned MaxRays)
 165 {
 166     *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 167     *poly = e = &e->x.p->arr[2*i+1];
 168     if (value_notzero_p(e->d))
 169         return;
 170     if (e->x.p->type != relation)
 171         return;
 172     if (e->x.p->size > 2)
 173         return;
 174     evalue *fr = &e->x.p->arr[0];
 175     assert(value_zero_p(fr->d));
 176     assert(fr->x.p->type == fractional);
 177     assert(fr->x.p->size == 3);
 178     Matrix *T = Matrix_Alloc(2, (*D)->Dimension+1);
 179     value_set_si(T->p[1][(*D)->Dimension], 1);
 180
 181     /* convert argument of fractional to polylib */
 182     /* the argument is assumed to be linear */
 183     evalue *p = &fr->x.p->arr[0];
 184     evalue_denom(p, &T->p[1][(*D)->Dimension]);
 185     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
 186         assert(p->x.p->type == polynomial);
 187         assert(p->x.p->size == 2);
 188         assert(value_notzero_p(p->x.p->arr[1].d));
 189         int pos = p->x.p->pos - 1;
 190         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
 191         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 192         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
 193     }
 194     int pos = (*D)->Dimension;
 195     value_assign(T->p[0][pos], p->x.n);
 196     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 197     value_division(T->p[0][pos], T->p[0][pos], p->d);
 198
 199     Polyhedron *E = NULL;
 200     for (Polyhedron *P = *D; P; P = P->next) {
 201         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
 202         I = DomainConstraintSimplify(I, MaxRays);
 203         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
 204         Polyhedron_Free(I);
 205         Polyhedron *next = P->next;
 206         P->next = NULL;
 207         Polyhedron *S = DomainIntersection(P, R, MaxRays);
 208         Polyhedron_Free(R);
 209         P->next = next;
 210         if (emptyQ2(S))
 211             Polyhedron_Free(S);
 212         else
 213             E = DomainConcat(S, E);
 214     }
 215     Matrix_Free(T);
 216
 217     *D = E;
 218     *poly = &e->x.p->arr[1];
 219 }
 220
 221 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 222                                       Polyhedron *ctx, const exvector& params)
 223 {
 224     piecewise_lst *pl;
 225     barvinok_options *options = barvinok_options_new_with_defaults();
 226     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
 227     barvinok_options_free(options);
 228     return pl;
 229 }
 230
 231 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 232                                       Polyhedron *ctx, const exvector& params,
 233                                       barvinok_options *options)
 234 {
 235     unsigned nparam = ctx->Dimension;
 236     if (EVALUE_IS_ZERO(*e))
 237         return pl_all;
 238     assert(value_zero_p(e->d));
 239     assert(e->x.p->type == partition);
 240     assert(e->x.p->size >= 2);
 241     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 242     unsigned PP_MaxRays = options->MaxRays;
 243     if (PP_MaxRays & POL_NO_DUAL)
 244         PP_MaxRays = 0;
 245
 246     exvector vars = constructVariableVector(nvars, "v");
 247     exvector allvars = vars;
 248     allvars.insert(allvars.end(), params.begin(), params.end());
 249
 250     for (int i = 0; i < e->x.p->size/2; ++i) {
 251         Param_Polyhedron *PP;
 252         Polyhedron *E;
 253         evalue *EP;
 254         Matrix *M;
 255         exvector floorvar;
 256
 257         evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
 258         ex poly = evalue2ex(EP, allvars, floorvar, &M);
 259         floorvar.insert(floorvar.end(), vars.begin(), vars.end());
 260         if (M) {
 261             Polyhedron *AE = align_context(E, M->NbColumns-2, options->MaxRays);
 262             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 263                 Domain_Free(E);
 264             E = DomainAddConstraints(AE, M, options->MaxRays);
 265             Matrix_Free(M);
 266             Domain_Free(AE);
 267         }
 268         if (is_exactly_a<fail>(poly)) {
 269             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 270                 Domain_Free(E);
 271             delete pl_all;
 272             return NULL;
 273         }
 274         for (Polyhedron *P = E; P; P = P->next) {
 275             Polyhedron *next = P->next;
 276             piecewise_lst *pl = new piecewise_lst(params);
 277             Polyhedron *P1 = P;
 278             P->next = NULL;
 279             PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 280             for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 281                 matrix VM = domainVertices(PP, Q, params);
 282                 lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 283                 pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
 284             }
 285             Param_Polyhedron_Free(PP);
 286             if (!pl_all)
 287                 pl_all = pl;
 288             else {
 289                 pl_all->combine(*pl);
 290                 delete pl;
 291             }
 292             P->next = next;
 293         }
 294         if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 295             Domain_Free(E);
 296     }
 297     return pl_all;
 298 }
 299
 300 }