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                       Vector *coset)
  69 {
  70     if (value_notzero_p(e->d))
  71         return value2numeric(e->x.n)/value2numeric(e->d);
  72     ex base_var = 0;
  73     ex poly = 0;
  74
  75     switch (e->x.p->type) {
  76     case polynomial:
  77         base_var = vars[e->x.p->pos-1];
  78         break;
  79     case flooring:
  80         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, false);
  81         break;
  82     case fractional:
  83         base_var = evalue2ex_add_var(&e->x.p->arr[0], extravar, expr, true);
  84         break;
  85     case periodic:
  86         assert(coset);
  87         return evalue2ex_r(&e->x.p->arr[VALUE_TO_INT(coset->p[e->x.p->pos-1])],
  88                            vars, extravar, expr, coset);
  89     default:
  90         return fail();
  91     }
  92
  93     int offset = type_offset(e->x.p);
  94     for (int i = e->x.p->size-1; i >= offset; --i) {
  95         poly *= base_var;
  96         ex t = evalue2ex_r(&e->x.p->arr[i], vars, extravar, expr, coset);
  97         if (is_exactly_a<fail>(t))
  98             return t;
  99         poly += t;
 100     }
 101     return poly;
 102 }
 103
 104 /* For each t = floor(e/d), set up two constraints
 105  *
 106  *       e - d t       >= 0
 107  *      -e + d t + d-1 >= 0
 108  *
 109  * e is assumed to be an affine expression.
 110  *
 111  * For each t = fract(e/d), set up two constraints
 112  *
 113  *          -d t + d-1 >= 0
 114  *             t       >= 0
 115  */
 116 static Matrix *setup_constraints(const vector<typed_evalue> expr, int nvar)
 117 {
 118     int extra = expr.size();
 119     if (!extra)
 120         return NULL;
 121     Matrix *M = Matrix_Alloc(2*extra, 1+extra+nvar+1);
 122     for (int i = 0; i < extra; ++i) {
 123         Value *d = &M->p[2*i][1+i];
 124         value_set_si(*d, 1);
 125         evalue_denom(expr[i].second, d);
 126         if (expr[i].first) {
 127             value_set_si(M->p[2*i][0], 1);
 128             value_decrement(M->p[2*i][1+extra+nvar], *d);
 129             value_oppose(*d, *d);
 130             value_set_si(M->p[2*i+1][0], 1);
 131             value_set_si(M->p[2*i+1][1+i], 1);
 132         } else {
 133             const evalue *e;
 134             for (e = expr[i].second; value_zero_p(e->d); e = &e->x.p->arr[0]) {
 135                 assert(e->x.p->type == polynomial);
 136                 assert(e->x.p->size == 2);
 137                 evalue *c = &e->x.p->arr[1];
 138                 value_multiply(M->p[2*i][1+extra+e->x.p->pos-1], *d, c->x.n);
 139                 value_division(M->p[2*i][1+extra+e->x.p->pos-1],
 140                                M->p[2*i][1+extra+e->x.p->pos-1], c->d);
 141             }
 142             value_multiply(M->p[2*i][1+extra+nvar], *d, e->x.n);
 143             value_division(M->p[2*i][1+extra+nvar], M->p[2*i][1+extra+nvar], e->d);
 144             value_oppose(*d, *d);
 145             value_set_si(M->p[2*i][0], -1);
 146             Vector_Scale(M->p[2*i], M->p[2*i+1], M->p[2*i][0], 1+extra+nvar+1);
 147             value_set_si(M->p[2*i][0], 1);
 148             value_subtract(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar], *d);
 149             value_decrement(M->p[2*i+1][1+extra+nvar], M->p[2*i+1][1+extra+nvar]);
 150         }
 151     }
 152     return M;
 153 }
 154
 155 static bool evalue_is_periodic(const evalue *e, Vector *periods)
 156 {
 157     int i, offset;
 158     bool is_periodic = false;
 159
 160     if (value_notzero_p(e->d))
 161         return false;
 162
 163     assert(e->x.p->type != partition);
 164     if (e->x.p->type == periodic) {
 165         Value size;
 166         value_init(size);
 167         value_set_si(size, e->x.p->size);
 168         value_lcm(periods->p[e->x.p->pos-1], size, &periods->p[e->x.p->pos-1]);
 169         value_clear(size);
 170         is_periodic = true;
 171     }
 172     offset = type_offset(e->x.p);
 173     for (i = e->x.p->size-1; i >= offset; --i)
 174         is_periodic = evalue_is_periodic(&e->x.p->arr[i], periods) || is_periodic;
 175     return is_periodic;
 176 }
 177
 178 static ex evalue2lst(const evalue *e, const exvector& vars,
 179                      exvector& extravar, vector<typed_evalue>& expr,
 180                      Vector *periods)
 181 {
 182     Vector *coset = Vector_Alloc(periods->Size);
 183     lst list;
 184     for (;;) {
 185         int i;
 186         list.append(evalue2ex_r(e, vars, extravar, expr, coset));
 187         for (i = coset->Size-1; i >= 0; --i) {
 188             value_increment(coset->p[i], coset->p[i]);
 189             if (value_lt(coset->p[i], periods->p[i]))
 190                 break;
 191             value_set_si(coset->p[i], 0);
 192         }
 193         if (i < 0)
 194             break;
 195     }
 196     Vector_Free(coset);
 197     return list;
 198 }
 199
 200 ex evalue2ex(const evalue *e, const exvector& vars, exvector& floorvar,
 201              Matrix **C, Vector **p)
 202 {
 203     vector<typed_evalue> expr;
 204     Vector *periods = Vector_Alloc(vars.size());
 205     assert(p);
 206     assert(C);
 207     for (int i = 0; i < periods->Size; ++i)
 208         value_set_si(periods->p[i], 1);
 209     if (evalue_is_periodic(e, periods)) {
 210         *p = periods;
 211         *C = NULL;
 212         lst list;
 213         return list;
 214     } else {
 215         Vector_Free(periods);
 216         *p = NULL;
 217         ex poly = evalue2ex_r(e, vars, floorvar, expr, NULL);
 218         Matrix *M = setup_constraints(expr, vars.size());
 219         *C = M;
 220         return poly;
 221     }
 222 }
 223
 224 /* if the evalue is a relation, we use the relation to cut off the
 225  * the edges of the domain
 226  */
 227 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
 228                                 unsigned MaxRays)
 229 {
 230     *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 231     *poly = e = &e->x.p->arr[2*i+1];
 232     if (value_notzero_p(e->d))
 233         return;
 234     if (e->x.p->type != relation)
 235         return;
 236     if (e->x.p->size > 2)
 237         return;
 238     evalue *fr = &e->x.p->arr[0];
 239     assert(value_zero_p(fr->d));
 240     assert(fr->x.p->type == fractional);
 241     assert(fr->x.p->size == 3);
 242     Matrix *T = Matrix_Alloc(2, (*D)->Dimension+1);
 243     value_set_si(T->p[1][(*D)->Dimension], 1);
 244
 245     /* convert argument of fractional to polylib */
 246     /* the argument is assumed to be linear */
 247     evalue *p = &fr->x.p->arr[0];
 248     evalue_denom(p, &T->p[1][(*D)->Dimension]);
 249     for (;value_zero_p(p->d); p = &p->x.p->arr[0]) {
 250         assert(p->x.p->type == polynomial);
 251         assert(p->x.p->size == 2);
 252         assert(value_notzero_p(p->x.p->arr[1].d));
 253         int pos = p->x.p->pos - 1;
 254         value_assign(T->p[0][pos], p->x.p->arr[1].x.n);
 255         value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 256         value_division(T->p[0][pos], T->p[0][pos], p->x.p->arr[1].d);
 257     }
 258     int pos = (*D)->Dimension;
 259     value_assign(T->p[0][pos], p->x.n);
 260     value_multiply(T->p[0][pos], T->p[0][pos], T->p[1][(*D)->Dimension]);
 261     value_division(T->p[0][pos], T->p[0][pos], p->d);
 262
 263     Polyhedron *E = NULL;
 264     for (Polyhedron *P = *D; P; P = P->next) {
 265         Polyhedron *I = Polyhedron_Image(P, T, MaxRays);
 266         I = DomainConstraintSimplify(I, MaxRays);
 267         Polyhedron *R = Polyhedron_Preimage(I, T, MaxRays);
 268         Polyhedron_Free(I);
 269         Polyhedron *next = P->next;
 270         P->next = NULL;
 271         Polyhedron *S = DomainIntersection(P, R, MaxRays);
 272         Polyhedron_Free(R);
 273         P->next = next;
 274         if (emptyQ2(S))
 275             Polyhedron_Free(S);
 276         else
 277             E = DomainConcat(S, E);
 278     }
 279     Matrix_Free(T);
 280
 281     *D = E;
 282     *poly = &e->x.p->arr[1];
 283 }
 284
 285 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 286                                       Polyhedron *ctx, const exvector& params)
 287 {
 288     piecewise_lst *pl;
 289     barvinok_options *options = barvinok_options_new_with_defaults();
 290     pl = evalue_bernstein_coefficients(pl_all, e, ctx, params, options);
 291     barvinok_options_free(options);
 292     return pl;
 293 }
 294
 295 static piecewise_lst *bernstein_coefficients(piecewise_lst *pl_all,
 296                             Polyhedron *D, const ex& poly,
 297                             Polyhedron *ctx,
 298                             const exvector& params, const exvector& floorvar,
 299                             barvinok_options *options)
 300 {
 301     if (emptyQ2(D))
 302         return pl_all;
 303
 304     unsigned PP_MaxRays = options->MaxRays;
 305     if (PP_MaxRays & POL_NO_DUAL)
 306         PP_MaxRays = 0;
 307
 308     for (Polyhedron *P = D; P; P = P->next) {
 309         Param_Polyhedron *PP;
 310         Polyhedron *next = P->next;
 311         Polyhedron *P1 = P;
 312         P->next = NULL;
 313         if (Polyhedron_is_unbounded(P, ctx->Dimension, options->MaxRays)) {
 314             fprintf(stderr, "warning: unbounded domain skipped\n");
 315             Polyhedron_Print(stderr, P_VALUE_FMT, P);
 316             P->next = next;
 317             continue;
 318         }
 319         PP = Polyhedron2Param_Domain(P, ctx, PP_MaxRays);
 320         piecewise_lst *pl = new piecewise_lst(params);
 321         for (Param_Domain *Q = PP->D; Q; Q = Q->next) {
 322             matrix VM = domainVertices(PP, Q, params);
 323             lst coeffs = bernsteinExpansion(VM, poly, floorvar, params);
 324             pl->list.push_back(guarded_lst(Polyhedron_Copy(Q->Domain), coeffs));
 325         }
 326         Param_Polyhedron_Free(PP);
 327         if (!pl_all)
 328             pl_all = pl;
 329         else {
 330             pl_all->combine(*pl);
 331             delete pl;
 332         }
 333         P->next = next;
 334     }
 335     return pl_all;
 336 }
 337
 338 /* Compute the coefficients of the polynomial corresponding to each coset
 339  * on its own domain.  This allows us to cut the domain on multiples of
 340  * the period.
 341  * To perform the cutting for a coset "i mod n = c" we map the domain
 342  * to the quotient space trough "i = i' n + c", simplify the constraints
 343  * (implicitly) and then map back to the original space.
 344  */
 345 static piecewise_lst *bernstein_coefficients_periodic(piecewise_lst *pl_all,
 346                             Polyhedron *D, const evalue *e,
 347                             Polyhedron *ctx, const exvector& vars,
 348                             const exvector& params, Vector *periods,
 349                             barvinok_options *options)
 350 {
 351     assert(D->Dimension == periods->Size);
 352     Matrix *T = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 353     Matrix *T2 = Matrix_Alloc(D->Dimension+1, D->Dimension+1);
 354     Vector *coset = Vector_Alloc(periods->Size);
 355     exvector extravar;
 356     vector<typed_evalue> expr;
 357     exvector allvars = vars;
 358     allvars.insert(allvars.end(), params.begin(), params.end());
 359
 360     value_set_si(T2->p[D->Dimension][D->Dimension], 1);
 361     for (int i = 0; i < D->Dimension; ++i) {
 362         value_assign(T->p[i][i], periods->p[i]);
 363         value_lcm(T2->p[D->Dimension][D->Dimension], periods->p[i],
 364                   &T2->p[D->Dimension][D->Dimension]);
 365     }
 366     value_set_si(T->p[D->Dimension][D->Dimension], 1);
 367     for (int i = 0; i < D->Dimension; ++i)
 368         value_division(T2->p[i][i], T2->p[D->Dimension][D->Dimension],
 369                        periods->p[i]);
 370     for (;;) {
 371         int i;
 372         ex poly = evalue2ex_r(e, allvars, extravar, expr, coset);
 373         assert(extravar.size() == 0);
 374         assert(expr.size() == 0);
 375         Polyhedron *E = DomainPreimage(D, T, options->MaxRays);
 376         Polyhedron *F = DomainPreimage(E, T2, options->MaxRays);
 377         Polyhedron_Free(E);
 378         if (!emptyQ2(F))
 379             pl_all = bernstein_coefficients(pl_all, F, poly, ctx, params,
 380                                             vars, options);
 381         Polyhedron_Free(F);
 382         for (i = D->Dimension-1; i >= 0; --i) {
 383             value_increment(coset->p[i], coset->p[i]);
 384             value_increment(T->p[i][D->Dimension], T->p[i][D->Dimension]);
 385             value_subtract(T2->p[i][D->Dimension], T2->p[i][D->Dimension],
 386                            T2->p[i][i]);
 387             if (value_lt(coset->p[i], periods->p[i]))
 388                 break;
 389             value_set_si(coset->p[i], 0);
 390             value_set_si(T->p[i][D->Dimension], 0);
 391             value_set_si(T2->p[i][D->Dimension], 0);
 392         }
 393         if (i < 0)
 394             break;
 395     }
 396     Vector_Free(coset);
 397     Matrix_Free(T);
 398     Matrix_Free(T2);
 399     return pl_all;
 400 }
 401
 402 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
 403                                       Polyhedron *ctx, const exvector& params,
 404                                       barvinok_options *options)
 405 {
 406     unsigned nparam = ctx->Dimension;
 407     if (EVALUE_IS_ZERO(*e))
 408         return pl_all;
 409     assert(value_zero_p(e->d));
 410     assert(e->x.p->type == partition);
 411     assert(e->x.p->size >= 2);
 412     unsigned nvars = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension - nparam;
 413
 414     exvector vars = constructVariableVector(nvars, "v");
 415     exvector allvars = vars;
 416     allvars.insert(allvars.end(), params.begin(), params.end());
 417
 418     for (int i = 0; i < e->x.p->size/2; ++i) {
 419         Polyhedron *E;
 420         evalue *EP;
 421         Matrix *M;
 422         Vector *periods;
 423         exvector floorvar;
 424
 425         evalue_extract_poly(e, i, &E, &EP, options->MaxRays);
 426         ex poly = evalue2ex(EP, allvars, floorvar, &M, &periods);
 427         floorvar.insert(floorvar.end(), vars.begin(), vars.end());
 428         if (M) {
 429             Polyhedron *AE = align_context(E, M->NbColumns-2, options->MaxRays);
 430             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 431                 Domain_Free(E);
 432             E = DomainAddConstraints(AE, M, options->MaxRays);
 433             Matrix_Free(M);
 434             Domain_Free(AE);
 435         }
 436         if (is_exactly_a<fail>(poly)) {
 437             if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 438                 Domain_Free(E);
 439             delete pl_all;
 440             return NULL;
 441         }
 442         if (periods)
 443             pl_all = bernstein_coefficients_periodic(pl_all, E, EP, ctx, vars,
 444                                                      params, periods, options);
 445         else
 446             pl_all = bernstein_coefficients(pl_all, E, poly, ctx, params,
 447                                             floorvar, options);
 448         if (periods)
 449             Vector_Free(periods);
 450         if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
 451             Domain_Free(E);
 452     }
 453     return pl_all;
 454 }
 455
 456 }