12 #include <bernstein/bernstein.h>
13 #include "polynomial.h"
16 using namespace GiNaC
;
22 static std::string
int2String(int n
);
23 static unsigned int findMaxDegree(ex polynomial
, const exvector
& Vars
);
24 static matrix
getAiMatrix(unsigned int nbVert
);
25 static ex
getBasis(unsigned int nbVert
, matrix
&A
);
26 static ex
powerMonomials(polynomial
&poly
, matrix
&A
, unsigned int nbVert
,
27 unsigned int maxDegree
, ex
&basis
);
28 static lst
getCoefficients(ex hom
, unsigned d
, const matrix
& A
);
30 exvector
constructParameterVector(const char * const *param_names
, unsigned nbParams
)
33 for (int i
= 0; i
< nbParams
; ++i
) {
34 P
[i
] = symbol(param_names
[i
]);
36 cout
<< "P: " << P
[i
] << endl
;
43 exvector
constructVariableVector(unsigned nbVariables
, const char *prefix
)
45 exvector
V(nbVariables
);
46 for (int i
= 0; i
< nbVariables
; ++i
) {
47 V
[i
] = symbol(prefix
+ int2String(i
));
49 cout
<< "V: " << V
[i
] << endl
;
56 numeric
value2numeric(const Value v
)
58 int sa
= v
[0]._mp_size
;
61 int abs_sa
= sa
< 0 ? -sa
: sa
;
63 for (int i
= abs_sa
-1; i
>= 0; --i
) {
64 res
= res
<< GMP_LIMB_BITS
;
65 res
= res
+ v
[0]._mp_d
[i
];
67 return numeric(sa
< 0 ? -res
: res
);
71 #if (GMP_LIMB_BITS==32)
72 inline mp_limb_t
cl_I_to_limb(const cln::cl_I
& x
) { return cln::cl_I_to_UL(x
); }
73 #elif (GMP_LIMB_BITS==64)
74 inline mp_limb_t
cl_I_to_limb(const cln::cl_I
& x
) { return cln::cl_I_to_UQ(x
); }
77 void numeric2value(const numeric
& n
, Value
& v
)
80 cln::cl_I abs_n
= cln::the
<cln::cl_I
>(abs(n
).to_cl_N());
83 for (abs_sa
= 0; abs_n
!= 0; abs_sa
++)
84 abs_n
= abs_n
>> GMP_LIMB_BITS
;
85 _mpz_realloc(v
, abs_sa
);
88 mask
= mask
<< GMP_LIMB_BITS
;
90 abs_n
= cln::the
<cln::cl_I
>(abs(n
).to_cl_N());
91 for (int i
= 0; i
< abs_sa
; ++i
) {
92 cln::cl_I digit
= abs_n
& mask
;
93 v
[0]._mp_d
[i
] = cl_I_to_limb(digit
);
94 abs_n
>> GMP_LIMB_BITS
;
97 v
[0]._mp_size
= n
< 0 ? -abs_sa
: abs_sa
;
101 matrix
domainVertices(Param_Polyhedron
*PP
, Param_Domain
*Q
, const exvector
& params
)
104 unsigned nbVertices
= 0;
105 unsigned nbRows
, nbColumns
;
109 nbRows
= PP
->V
->Vertex
->NbRows
;
110 nbColumns
= PP
->V
->Vertex
->NbColumns
;
112 FORALL_PVertex_in_ParamPolyhedron(V
, Q
, PP
)
114 END_FORALL_PVertex_in_ParamPolyhedron
;
116 matrix
VM(nbVertices
, nbRows
); // vertices matrix
119 FORALL_PVertex_in_ParamPolyhedron(V
, Q
, PP
)
120 for (unsigned i
= 0; i
< nbRows
; i
++) {
122 for (unsigned j
= 0; j
< nbColumns
-2; j
++)
123 t
+= value2numeric(V
->Vertex
->p
[i
][j
]) * params
[j
];
124 t
+= value2numeric(V
->Vertex
->p
[i
][nbColumns
-2]);
125 t
/= value2numeric(V
->Vertex
->p
[i
][nbColumns
-1]);
127 cout
<< "T: " << t
<< endl
;
132 END_FORALL_PVertex_in_ParamPolyhedron
;
138 * Do the Bernstein Expansion
140 * vert: vertices matrix
141 * poly: polynomial expression
142 * vars: vector of variables
143 * Params: vector of parameter
145 lst
bernsteinExpansion(const matrix
& vert
, const ex
& poly
, const exvector
& vars
,
146 const exvector
& Params
)
148 if (is_exactly_a
<lst
>(poly
)) {
149 lst polylst
= ex_to
<lst
>(poly
);
150 lst::const_iterator j
= polylst
.begin();
152 lst coeff
= bernsteinExpansion(vert
, *j
, vars
, Params
);
154 for (++j
; j
!= polylst
.end(); ++j
) {
155 lst::const_iterator k
;
156 lst new_coeff
= bernsteinExpansion(vert
, *j
, vars
, Params
);
157 for (k
= new_coeff
.begin(); k
!= new_coeff
.end(); ++k
)
159 coeff
.sort().unique();
164 unsigned maxDegree
= findMaxDegree(poly
, vars
);
165 matrix A
= getAiMatrix(vert
.rows());
168 cout
<< endl
<< "Polynomial: " << poly
<< endl
<< endl
;
169 cout
<< vert
<< endl
<< endl
;
170 cout
<< A
<< endl
<< endl
;
173 // obtain the variables value on the basis and replace
174 ex substitutions
= evalm(A
* vert
);
175 ex polynom
= replaceVariablesInPolynomial(poly
, vars
, substitutions
);
178 cout
<< variables
<< endl
<< endl
;
179 cout
<< "Preliminar Expansion: " << polynom
<< endl
<< endl
;
182 ex basis
= getBasis(vert
.rows(), A
);
185 cout
<< "Basis: " << basis
<< endl
<< endl
;
188 // monomials to n degree
189 polynomial
p(polynom
);
190 ex maxDegreePolynomial
= powerMonomials(p
, A
, vert
.rows(), maxDegree
, basis
);
192 return getCoefficients(maxDegreePolynomial
, maxDegree
, A
);
196 * Construct A_i matrix
198 * nbVert: number of vertices
200 matrix
getAiMatrix(unsigned int nbVert
)
202 matrix
A(1, nbVert
); // a_i matrix
203 for(unsigned int i
= 0; i
< nbVert
; i
++) {
204 A(0,i
) = symbol("a" + int2String(i
));
206 cout
<< "A: " << A(0,i
) << endl
;
214 * Construct the basis
217 * nbVert: number of vertices
219 ex
getBasis(unsigned int nbVert
, matrix
&A
)
222 for(unsigned int i
= 0; i
< nbVert
; i
++) {
230 * Finds the less than maxDegree monomials and multiply them
233 * polynomial: original polynomial
235 * nbVert: number of vertices
236 * maxDegree: maximal polynomial multidegree
237 * basis: basis of the polytope
239 ex
powerMonomials(polynomial
&poly
, matrix
&A
, unsigned int nbVert
240 , unsigned int maxDegree
, ex
&basis
)
242 ex maxDegreePolynomial
;
244 cout
<< "- Degree --------------------------------------" << endl
;
246 for (size_t i
= 0; i
!= poly
.nbTerms(); ++i
) {
247 unsigned int degree
= 0; // degree of the monomial
249 for(unsigned int j
= 0; j
< nbVert
; j
++) {
250 degree
+= poly
.term(i
).degree(A(0,j
));
253 cout
<< poly
.term(i
) << " Degree: " << degree
;
255 if(degree
< maxDegree
) {
256 ex degreeUp
= poly
.term(i
) * pow(basis
, maxDegree
- degree
);
258 cout
<< " --> New Term: " << degreeUp
.expand();
260 maxDegreePolynomial
+= degreeUp
.expand();
262 maxDegreePolynomial
+= poly
.term(i
);
265 cout
<< endl
<< "-----------------------------------------------" << endl
;
270 cout
<< endl
<< "Final Expansion: " << maxDegreePolynomial
<< endl
<< endl
;
272 return maxDegreePolynomial
;
278 * Finds the coefficients of the polynomial in terms of the Bernstein basis
280 * hom: homogeneous polynomial of degree d
283 * For each monomial of multi-degree (k[0], k[1], ..., k[n-1])
284 * we divide the corresponding coefficient by the multinomial
285 * coefficient d!/k[0]! k[1]! ... k[n-1]!
287 * The code looks a bit complicated because it's an iterative
288 * implementation of a recursive procedure.
289 * For each variable from 0 to n-1, we loop over the possible
290 * powers: 0..d for the first; 0..d-k[0]=left[0] for the second; ...
291 * Just for fun, we loop through these in the opposite direction
292 * only for the first variable.
294 * c[i] contains the coefficient of the selected powers of the first i+1 vars
295 * multinom[i] contains the partial multinomial coefficient.
297 lst
getCoefficients(ex hom
, unsigned d
, const matrix
& A
)
302 /* we should probably notice sooner that there is just one vertex */
304 coeff
.append(hom
.coeff(A(0, 0), d
));
315 for (k
[0] = d
; k
[0] >= 0; --k
[0]) {
316 c
[0] = hom
.coeff(A(0, 0), k
[0]);
320 multinom
[i
] = multinom
[0];
323 for (int j
= 2; j
<= left
[i
-1]; ++j
)
325 coeff
.append(c
[i
-1].coeff(A(0, i
), left
[i
-1]) /
330 if (k
[i
] >= left
[i
-1]) {
337 c
[i
] = c
[i
-1].coeff(A(0, i
), k
[i
]);
338 left
[i
] = left
[i
-1] - k
[i
];
340 multinom
[i
+1] = multinom
[i
];
345 return coeff
.sort().unique();
349 // replace the variables in the polynomial
350 ex
replaceVariablesInPolynomial(const ex
&poly
, const exvector
&vars
,
351 const ex
&substitutions
)
355 for(unsigned int i
= 0; i
< vars
.size(); i
++) {
357 cout
<< "Replacing: " << vars
[i
] << " by " << substitutions
[i
] << endl
;
359 replace
.append(vars
[i
] == substitutions
[i
]);
361 ex polyRepl
= poly
.subs(replace
);
363 return(polyRepl
.expand());
367 /* Converts int n to string */
368 string
int2String(int n
)
370 char numeroVariable
[DIGITS
];
371 snprintf(numeroVariable
, DIGITS
, "%d", n
);
372 string
nroV(numeroVariable
);
379 * Find the maximum multi-degree of the polinomial
381 * polynomial: polynomial
382 * Vars: variables matrix
383 * nbVar: number of variables
385 static unsigned int findMaxDegree(ex polynomial
, const exvector
& Vars
, int pos
)
387 if (pos
>= Vars
.size())
391 for (int i
= 0; i
<= polynomial
.degree(Vars
[pos
]); ++i
) {
392 unsigned degree
= i
+ findMaxDegree(polynomial
.coeff(Vars
[pos
], i
),
400 unsigned int findMaxDegree(ex polynomial
, const exvector
& Vars
)
402 return findMaxDegree(polynomial
, Vars
, 0);