barvinok_maximize: new tool for maximizing piecewise quasi-polynomial
[barvinok.git] / bernstein.cc
blob0f726999669864eed72a09f2aec764f06f91a68a
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>
8 using namespace GiNaC;
9 using namespace bernstein;
11 namespace barvinok {
13 ex evalue2ex(evalue *e, const exvector& vars)
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;
27 return poly;
30 /* if the evalue is a relation, we use the relation to cut off the
31 * the edges of the domain
33 static void evalue_extract_poly(evalue *e, int i, Polyhedron **D, evalue **poly,
34 unsigned MaxRays)
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);
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);
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);
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);
85 Matrix_Free(T);
87 *D = E;
88 *poly = &e->x.p->arr[1];
91 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
92 Polyhedron *ctx, const exvector& params)
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;
101 piecewise_lst *evalue_bernstein_coefficients(piecewise_lst *pl_all, evalue *e,
102 Polyhedron *ctx, const exvector& params,
103 barvinok_options *options)
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;
116 exvector vars = constructVariableVector(nvars, "v");
117 exvector allvars = vars;
118 allvars.insert(allvars.end(), params.begin(), params.end());
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;
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));
143 Param_Polyhedron_Free(PP);
144 if (!pl_all)
145 pl_all = pl;
146 else {
147 pl_all->combine(*pl);
148 delete pl;
150 P->next = next;
152 if (E != EVALUE_DOMAIN(e->x.p->arr[2*i]))
153 Domain_Free(E);
155 return pl_all;