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