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>
9 using namespace bernstein
;
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
)
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
))
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
,
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
))
40 if (e
->x
.p
->type
!= relation
)
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
);
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
);
75 Polyhedron
*next
= P
->next
;
77 Polyhedron
*S
= DomainIntersection(P
, R
, MaxRays
);
83 E
= DomainConcat(S
, 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
)
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
);
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
))
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
)
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
;
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
]))
132 for (Polyhedron
*P
= E
; P
; P
= P
->next
) {
133 Polyhedron
*next
= P
->next
;
134 piecewise_lst
*pl
= new piecewise_lst(params
);
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
);
147 pl_all
->combine(*pl
);
152 if (E
!= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]))