bernstein/src/bernstein.cpp

   1 /*
   2  *      Bernstein Expansion
   3  *
   4  *      - ginac functions
   5  */
   6
   7 #include <assert.h>
   8 #include <iostream>
   9 #include <cln/cln.h>
  10 #include <string>
  11
  12 #include <bernstein/bernstein.h>
  13 #include "polynomial.h"
  14
  15 using namespace std;
  16 using namespace GiNaC;
  17
  18 #define DIGITS 10
  19
  20 namespace bernstein {
  21
  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);
  29
  30 exvector constructParameterVector(const char * const *param_names, unsigned nbParams)
  31 {
  32         exvector P(nbParams);
  33         for (int i = 0; i < nbParams; ++i) {
  34                 P[i] = symbol(param_names[i]);
  35 #ifdef DEBUG
  36                 cout << "P: " << P[i] << endl;
  37 #endif
  38         }
  39         return P;
  40 }
  41
  42
  43 exvector constructVariableVector(unsigned nbVariables, const char *prefix)
  44 {
  45         exvector V(nbVariables);
  46         for (int i = 0; i < nbVariables; ++i) {
  47                 V[i] = symbol(prefix + int2String(i));
  48 #ifdef DEBUG
  49                 cout << "V: " << V[i] << endl;
  50 #endif
  51         }
  52         return V;
  53 }
  54
  55
  56 numeric value2numeric(const Value v)
  57 {
  58     int sa = v[0]._mp_size;
  59     if (!sa)
  60         return 0;
  61     int abs_sa = sa < 0 ? -sa : sa;
  62     cln::cl_I res = 0;
  63     for (int i = abs_sa-1; i >= 0; --i) {
  64         res = res << GMP_LIMB_BITS;
  65         res = res + v[0]._mp_d[i];
  66     }
  67     return numeric(sa < 0 ? -res : res);
  68 }
  69
  70
  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); }
  75 #endif
  76
  77 void numeric2value(const numeric& n, Value& v)
  78 {
  79     cln::cl_I mask;
  80     cln::cl_I abs_n = cln::the<cln::cl_I>(abs(n).to_cl_N());
  81     int abs_sa;
  82
  83     for (abs_sa = 0; abs_n != 0; abs_sa++)
  84         abs_n = abs_n >> GMP_LIMB_BITS;
  85     _mpz_realloc(v, abs_sa);
  86
  87     mask = 1;
  88     mask = mask << GMP_LIMB_BITS;
  89     mask = mask - 1;
  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 = abs_n >> GMP_LIMB_BITS;
  95     }
  96
  97     v[0]._mp_size = n < 0 ? -abs_sa : abs_sa;
  98 }
  99
 100
 101 matrix domainVertices(Param_Polyhedron *PP, Param_Domain *Q, const exvector& params)
 102 {
 103         Param_Vertices *V;
 104         unsigned nbVertices = 0;
 105         unsigned nbRows, nbColumns;
 106         int v;
 107
 108         assert(PP->nbV > 0);
 109         nbRows = PP->V->Vertex->NbRows;
 110         nbColumns = PP->V->Vertex->NbColumns;
 111
 112         FORALL_PVertex_in_ParamPolyhedron(V, Q, PP)
 113                 ++nbVertices;
 114         END_FORALL_PVertex_in_ParamPolyhedron;
 115
 116         matrix VM(nbVertices, nbRows);          // vertices matrix
 117
 118         v = 0;
 119         FORALL_PVertex_in_ParamPolyhedron(V, Q, PP)
 120                 for (unsigned i = 0; i < nbRows; i++) {
 121                         ex t;
 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]);
 126 #ifdef DEBUG
 127                         cout << "T: " << t << endl;
 128 #endif
 129                         VM(v, i) = t;
 130                 }
 131                 ++v;
 132         END_FORALL_PVertex_in_ParamPolyhedron;
 133
 134         return VM;
 135 }
 136
 137 /*
 138  * Do the Bernstein Expansion
 139  *
 140  *      vert: vertices matrix
 141  *      poly: polynomial expression
 142  *      vars: vector of variables
 143  *      Params: vector of parameter
 144  */
 145 lst bernsteinExpansion(const matrix& vert, const ex& poly, const exvector& vars,
 146                        const exvector& Params)
 147 {
 148         if (is_exactly_a<lst>(poly)) {
 149                 lst polylst = ex_to<lst>(poly);
 150                 lst::const_iterator j = polylst.begin();
 151
 152                 lst coeff = bernsteinExpansion(vert, *j, vars, Params);
 153
 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)
 158                             coeff.append(*k);
 159                         coeff.sort().unique();
 160                 }
 161                 return coeff;
 162         }
 163
 164         unsigned maxDegree = findMaxDegree(poly, vars);
 165         matrix A = getAiMatrix(vert.rows());
 166
 167 #ifdef DEBUG
 168         cout << endl << "Polynomial: " << poly << endl << endl;
 169         cout << vert << endl << endl;
 170         cout << A << endl << endl;
 171 #endif
 172
 173         // obtain the variables value on the basis and replace
 174         ex substitutions = evalm(A * vert);
 175         ex polynom = replaceVariablesInPolynomial(poly, vars, substitutions);
 176
 177 #ifdef DEBUG
 178         cout << substitutions << endl << endl;
 179         cout << "Preliminar Expansion: " << polynom << endl << endl;
 180 #endif
 181
 182         ex basis = getBasis(vert.rows(), A);
 183
 184 #ifdef DEBUG
 185         cout << "Basis: " << basis << endl<< endl;
 186 #endif
 187
 188         // monomials to n degree
 189         polynomial p(polynom);
 190         ex maxDegreePolynomial = powerMonomials(p, A, vert.rows(), maxDegree, basis);
 191
 192         return getCoefficients(maxDegreePolynomial, maxDegree, A);
 193 }
 194
 195 /*
 196  * Construct A_i matrix
 197  *
 198  *      nbVert: number of vertices
 199  */
 200 matrix getAiMatrix(unsigned int nbVert)
 201 {
 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));
 205 #ifdef DEBUG
 206                 cout << "A: " << A(0,i) << endl;
 207 #endif
 208         }
 209         return A;
 210 }
 211
 212
 213 /*
 214  * Construct the basis
 215  *
 216  *      A: a_i matrix
 217  *      nbVert: number of vertices
 218  */
 219 ex getBasis(unsigned int nbVert, matrix &A)
 220 {
 221         ex basis;
 222         for(unsigned int i = 0; i < nbVert; i++) {
 223                 basis += A(0,i);
 224         }
 225         return basis;
 226 }
 227
 228
 229 /*
 230  * Finds the less than maxDegree monomials and multiply them
 231  * by the basis
 232  *
 233  *      polynomial: original polynomial
 234  *      A: a_i matrix
 235  *      nbVert: number of vertices
 236  *      maxDegree: maximal polynomial multidegree
 237  *      basis: basis of the polytope
 238  */
 239 ex powerMonomials(polynomial &poly, matrix &A, unsigned int nbVert
 240                   , unsigned int maxDegree, ex &basis)
 241 {
 242         ex maxDegreePolynomial;
 243 #ifdef DEBUG
 244         cout << "- Degree --------------------------------------" << endl;
 245 #endif
 246         for (size_t i = 0; i != poly.nbTerms(); ++i) {
 247                 unsigned int degree = 0;                // degree of the monomial
 248
 249                 for(unsigned int j = 0; j < nbVert; j++) {
 250                         degree += poly.term(i).degree(A(0,j));
 251                 }
 252 #ifdef DEBUG
 253                 cout << poly.term(i) << " Degree: " << degree;
 254 #endif
 255                 if(degree < maxDegree) {
 256                         ex degreeUp = poly.term(i) * pow(basis, maxDegree - degree);
 257 #ifdef DEBUG
 258                         cout << "   --> New Term: " <<  degreeUp.expand();
 259 #endif
 260                         maxDegreePolynomial += degreeUp.expand();
 261                 } else {
 262                         maxDegreePolynomial += poly.term(i);
 263                 }
 264 #ifdef DEBUG
 265                 cout << endl << "-----------------------------------------------" << endl;
 266 #endif
 267
 268         }
 269 #ifdef DEBUG
 270         cout << endl << "Final Expansion: " << maxDegreePolynomial << endl << endl;
 271 #endif
 272         return maxDegreePolynomial;
 273 }
 274
 275
 276
 277 /*
 278  * Finds the coefficients of the polynomial in terms of the Bernstein basis
 279  *
 280  *      hom: homogeneous polynomial of degree d
 281  *      A: a_i matrix
 282  *
 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]!
 286  *
 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.
 293  *
 294  *      c[i] contains the coefficient of the selected powers of the first i+1 vars
 295  *      multinom[i] contains the partial multinomial coefficient.
 296  */
 297 lst getCoefficients(ex hom, unsigned d, const matrix& A)
 298 {
 299         lst coeff;
 300         int n = A.cols();
 301
 302         /* we should probably notice sooner that there is just one vertex */
 303         if (n == 1) {
 304             coeff.append(hom.coeff(A(0, 0), d));
 305             return coeff;
 306         }
 307
 308         ex c[n];
 309         int left[n];
 310         int k[n];
 311         numeric multinom[n];
 312         assert(n >= 2);
 313
 314         multinom[0] = 1;
 315         for (k[0] = d; k[0] >= 0; --k[0]) {
 316                 c[0] = hom.coeff(A(0, 0), k[0]);
 317                 left[0] = d - k[0];
 318                 int i = 1;
 319                 k[i] = -1;
 320                 multinom[i] = multinom[0];
 321                 while (i > 0) {
 322                         if (i == n-1) {
 323                                 for (int j = 2; j <= left[i-1]; ++j)
 324                                         multinom[i] /= j;
 325                                 coeff.append(c[i-1].coeff(A(0, i), left[i-1]) /
 326                                                 multinom[i]);
 327                                 --i;
 328                                 continue;
 329                         }
 330                         if (k[i] >= left[i-1]) {
 331                                 --i;
 332                                 continue;
 333                         }
 334                         ++k[i];
 335                         if (k[i])
 336                                 multinom[i] /= k[i];
 337                         c[i] = c[i-1].coeff(A(0, i), k[i]);
 338                         left[i] = left[i-1] - k[i];
 339                         k[i+1] = -1;
 340                         multinom[i+1] = multinom[i];
 341                         ++i;
 342                 }
 343                 multinom[0] *= k[0];
 344         }
 345         return coeff.sort().unique();
 346 }
 347
 348
 349 // replace the variables in the polynomial
 350 ex replaceVariablesInPolynomial(const ex &poly, const exvector &vars,
 351                                 const ex &substitutions)
 352 {
 353         lst replace;
 354
 355         for(unsigned int i = 0; i < vars.size(); i++) {
 356 #ifdef DEBUG
 357                 cout << "Replacing: " << vars[i] << " by " << substitutions[i] << endl;
 358 #endif
 359                 replace.append(vars[i] == substitutions[i]);
 360         }
 361         ex polyRepl = poly.subs(replace);
 362
 363         return(polyRepl.expand());
 364 }
 365
 366
 367 /* Converts int n to string */
 368 string int2String(int n)
 369 {
 370         char numeroVariable[DIGITS];
 371         snprintf(numeroVariable, DIGITS, "%d", n);
 372         string nroV(numeroVariable);
 373
 374         return nroV;
 375 }
 376
 377
 378 /*
 379  * Find the maximum multi-degree of the polinomial
 380  *
 381  *      polynomial: polynomial
 382  *      Vars: variables matrix
 383  *      nbVar: number of variables
 384  */
 385 static unsigned int findMaxDegree(ex polynomial, const exvector& Vars, int pos)
 386 {
 387         if (pos >= Vars.size())
 388                 return 0;
 389
 390         unsigned max = 0;
 391         for (int i = 0; i <= polynomial.degree(Vars[pos]); ++i) {
 392                 unsigned degree = i + findMaxDegree(polynomial.coeff(Vars[pos], i),
 393                                                     Vars, pos+1);
 394                 if (degree > max)
 395                         max = degree;
 396         }
 397         return max;
 398 }
 399
 400 unsigned int findMaxDegree(ex polynomial, const exvector& Vars)
 401 {
 402         return findMaxDegree(polynomial.expand(), Vars, 0);
 403 }
 404
 405 }