GetFEM  5.5
bgeot_poly.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2026 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 
32 #ifndef BGEOT_POLY_H__
33 #define BGEOT_POLY_H__
34 
35 /** @file bgeot_poly.h
36  @author Yves Renard <[email protected]>
37  @date December 01, 2000.
38  @brief Multivariate polynomials.
39 */
40 
41 #include "bgeot_config.h"
42 #include "dal_static_stored_objects.h"
43 #include <vector>
44 
45 namespace bgeot
46 {
47  /// used as the common size type in the library
48  typedef size_t size_type;
49  ///
50  /// used as the common short type integer in the library
51  typedef gmm::uint16_type short_type;
52  ///
53 
54  /** Return the value of @f$ \frac{(n+p)!}{n!p!} @f$ which
55  * is the number of monomials of a polynomial of @f$n@f$
56  * variables and degree @f$d@f$.
57  */
58  size_type alpha(short_type n, short_type d);
59 
60  /** Vector of integer (16 bits type) which represent the powers
61  * of a monomial
62  */
63  class power_index {
64  std::vector<short_type> v;
65  mutable short_type degree_;
66  mutable size_type global_index_;
67  void dirty() const
68  { degree_ = short_type(-1); global_index_ = size_type(-1); }
69  public :
70  typedef std::vector<short_type>::iterator iterator;
71  typedef std::vector<short_type>::const_iterator const_iterator;
72  typedef std::vector<short_type>::reverse_iterator reverse_iterator;
73  typedef std::vector<short_type>::const_reverse_iterator const_reverse_iterator;
74  short_type operator[](size_type idx) const { return v[idx]; }
75  short_type& operator[](size_type idx) { dirty(); return v[idx]; }
76 
77  iterator begin() { dirty(); return v.begin(); }
78  const_iterator begin() const { return v.begin(); }
79  iterator end() { dirty(); return v.end(); }
80  const_iterator end() const { return v.end(); }
81 
82  reverse_iterator rbegin() { dirty(); return v.rbegin(); }
83  const_reverse_iterator rbegin() const { return v.rbegin(); }
84  reverse_iterator rend() { dirty(); return v.rend(); }
85  const_reverse_iterator rend() const { return v.rend(); }
86 
87  size_type size() const { return v.size(); }
88  /// Gives the next power index
89  const power_index &operator ++();
90  /// Gives the next power index
91  const power_index operator ++(int)
92  { power_index res = *this; ++(*this); return res; }
93  /// Gives the next previous index
94  const power_index &operator --();
95  /// Gives the next previous index
96  const power_index operator --(int)
97  { power_index res = *this; --(*this); return res; }
98  /** Gives the global number of the index (i.e. the position of
99  * the corresponding monomial
100  */
101  size_type global_index() const;
102  /// Gives the degree.
103  short_type degree() const;
104  /// Constructor
105  power_index(short_type nn);
106  /// Constructor
107  power_index() { dirty(); }
108  };
109 
110  /**
111  * This class deals with plain polynomials with
112  * several variables.
113  *
114  * A polynomial of @f$n@f$ variables and degree @f$d@f$ is stored in a vector
115  * of @f$\alpha_d^n@f$ components.
116  *
117  * <h3>Example of code</h3>
118  *
119  * the following code is valid :
120  * @code
121  * #include<bgeot_poly.h>
122  * bgeot::polynomial<double> P, Q;
123  * P = bgeot::polynomial<double>(2,2,1); // P = x
124  * Q = bgeot::polynomial<double>(2,2,2); // Q = y
125  * P += Q; // P is equal to x+y.
126  * P *= Q; // P is equal to xy + y^2
127  * bgeot::power_index pi(P.dim());
128  * bgeot::polynomial<double>::const_iterator ite = Q.end();
129  * bgeot::polynomial<double>::const_iterator itb = Q.begin();
130  * for ( ; itb != ite; ++itb, ++pi)
131  * if (*itq != double(0))
132  * cout "there is x to the power " << pi[0]
133  * << " and y to the power "
134  * << pi[1] << " and a coefficient " << *itq << endl;
135  * @endcode
136  *
137  * <h3>Monomials ordering.</h3>
138  *
139  * The constant coefficient is placed first with the index 0.
140  * Two monomials of different degrees are ordered following
141  * there respective degree.
142  *
143  * If two monomials have the same degree, they are ordered with the
144  * degree of the mononomials without the n firsts variables which
145  * have the same degree. The index of the monomial
146  * @f$ x_0^{i_0}x_1^{i_1} ... x_{n-1}^{i_{n-1}} @f$
147  * is then
148  * @f$ \alpha_{d-1}^{n} + \alpha_{d-i_0-1}^{n-1}
149  * + \alpha_{d-i_0-i_1-1}^{n-2} + ... + \alpha_{i_{n-1}-1}^{1}, @f$
150  * where @f$d = \sum_{l=0}^{n-1} i_l@f$ is the degree of the monomial.
151  * (by convention @f$\alpha_{-1}^{n} = 0@f$).
152  *
153  * <h3>Dealing with the vector of power.</h3>
154  *
155  * The answer to the question : what is the next and previous
156  * monomial of @f$x_0^{i_0}x_1^{i_1} ... x_{n-1}^{i_{n-1}}@f$ in the
157  * vector is the following :
158  *
159  * To take the next coefficient, let @f$l@f$ be the last index between 0
160  * and @f$n-2@f$ such that @f$i_l \ne 0@f$ (@f$l = -1@f$ if there is not), then
161  * make the operations @f$a = i_{n-1}; i_{n-1} = 0; i_{l+1} = a+1;
162  * \mbox{ if } l \ge 0 \mbox{ then } i_l = i_l - 1@f$.
163  *
164  * To take the previous coefficient, let @f$l@f$ be the last index
165  * between 0 and @f$n-1@f$ such that @f$i_l \ne 0@f$ (if there is not, there
166  * is no previous monomial) then make the operations @f$a = i_l;
167  * i_l = 0; i_{n-1} = a - 1; \mbox{ if } l \ge 1 \mbox{ then }
168  * i_{l-1} = i_{l-1} + 1@f$.
169  *
170  * <h3>Direct product multiplication.</h3>
171  *
172  * This direct product multiplication of P and Q is the
173  * multiplication considering that the variables of Q follow the
174  * variables of P. The result is a polynomial with the number of
175  * variables of P plus the number of variables of Q.
176  * The resulting polynomials have a smaller degree.
177  *
178  */
179  template<typename T> class polynomial : public std::vector<T> {
180  protected :
181 
182  short_type n, d;
183 
184  public :
185 
186  typedef typename std::vector<T>::iterator iterator;
187  typedef typename std::vector<T>::const_iterator const_iterator;
188  typedef typename std::vector<T>::reverse_iterator reverse_iterator;
189  typedef typename std::vector<T>::const_reverse_iterator const_reverse_iterator;
190 
191  /// Gives the degree of the polynomial
192  short_type degree() const { return d; }
193  /** gives the degree of the polynomial, considering only non-zero
194  * coefficients
195  */
196  short_type real_degree() const;
197  /// Gives the dimension (number of variables)
198  short_type dim() const { return n; }
199  /// Change the degree of the polynomial to d.
200  void change_degree(short_type dd);
201  /** Add to the polynomial a monomial of coefficient a and
202  * correpsonding to the power index pi.
203  */
204  void add_monomial(const T &coeff, const power_index &power);
205  /// Add Q to P. P contains the result.
206  polynomial &operator +=(const polynomial &Q);
207  /// Subtract Q from P. P contains the result.
208  polynomial &operator -=(const polynomial &Q);
209  /// Add Q to P.
210  polynomial operator +(const polynomial &Q) const
211  { polynomial R = *this; R += Q; return R; }
212  /// Subtract Q from P.
213  polynomial operator -(const polynomial &Q) const
214  { polynomial R = *this; R -= Q; return R; }
215  polynomial operator -() const;
216  /// Multiply P with Q. P contains the result.
217  polynomial &operator *=(const polynomial &Q);
218  /// Multiply P with Q.
219  polynomial operator *(const polynomial &Q) const;
220  /** Product of P and Q considering that variables of Q come after
221  * variables of P. P contains the result
222  */
223  void direct_product(const polynomial &Q);
224  /// Multiply P with the scalar a. P contains the result.
225  polynomial &operator *=(const T &e);
226  /// Multiply P with the scalar a.
227  polynomial operator *(const T &e) const;
228  /// Divide P with the scalar a. P contains the result.
229  polynomial &operator /=(const T &e);
230  /// Divide P with the scalar a.
231  polynomial operator /(const T &e) const
232  { polynomial res = *this; res /= e; return res; }
233  /// operator ==.
234  bool operator ==(const polynomial &Q) const;
235  /// operator !=.
236  bool operator !=(const polynomial &Q) const
237  { return !(operator ==(*this,Q)); }
238  /// Derivative of P with respect to the variable k. P contains the result.
239  void derivative(short_type k);
240  /// Makes P = 1.
241  void one() { change_degree(0); (*this)[0] = T(1); }
242  void clear() { change_degree(0); (*this)[0] = T(0); }
243  bool is_zero()
244  { return(this->real_degree()==0) && (this->size()==0 || (*this)[0]==T(0)); }
245  template <typename ITER> T horner(power_index &mi, short_type k,
246  short_type de, const ITER &it) const;
247  /** Evaluate the polynomial. "it" is an iterator pointing to the list
248  * of variables. A Horner scheme is used.
249  */
250  template <typename ITER> T eval(const ITER &it) const;
251 
252  /// Constructor.
253  polynomial() : std::vector<T>(1)
254  { n = 0; d = 0; (*this)[0] = 0.0; }
255  /// Constructor.
256  polynomial(short_type dim_, short_type degree_);
257  /// Constructor for the polynomial 'x' (k=0), 'y' (k=1), 'z' (k=2) etc.
258  polynomial(short_type dim_, short_type degree_, short_type k);
259  };
260 
261 
262  template<typename T> polynomial<T>::polynomial(short_type nn, short_type dd)
263  : std::vector<T>(alpha(nn,dd))
264  { n = nn; d = dd; std::fill(this->begin(), this->end(), T(0)); }
265 
266  template<typename T> polynomial<T>::polynomial(short_type nn,
267  short_type dd, short_type k)
268  : std::vector<T>(alpha(nn,dd)) {
269  n = nn; d = std::max(short_type(1), dd);
270  std::fill(this->begin(), this->end(), T(0));
271  (*this)[k+1] = T(1);
272  }
273 
274  template<typename T>
275  polynomial<T> polynomial<T>::operator *(const polynomial &Q) const
276  { polynomial res = *this; res *= Q; return res; }
277 
278  template<typename T>
279  bool polynomial<T>::operator ==(const polynomial &Q) const {
280  if (dim() != Q.dim()) return false;
281  const_iterator it1 = this->begin(), ite1 = this->end();
282  const_iterator it2 = Q.begin(), ite2 = Q.end();
283  for ( ; it1 != ite1 && it2 != ite2; ++it1, ++it2)
284  if (*it1 != *it2) return false;
285  for ( ; it1 != ite1; ++it1) if (*it1 != T(0)) return false;
286  for ( ; it2 != ite2; ++it2) if (*it2 != T(0)) return false;
287  return true;
288  }
289 
290  template<typename T>
291  polynomial<T> polynomial<T>::operator *(const T &e) const
292  { polynomial res = *this; res *= e; return res; }
293 
294  template<typename T> short_type polynomial<T>::real_degree() const {
295  const_reverse_iterator it = this->rbegin(), ite = this->rend();
296  size_type l = this->size();
297  for ( ; it != ite; ++it, --l) { if (*it != T(0)) break; }
298  short_type dd = degree();
299  while (dd > 0 && alpha(n, short_type(dd-1)) > l) --dd;
300  return dd;
301  }
302 
303  template<typename T> void polynomial<T>::change_degree(short_type dd) {
304  this->resize(alpha(n,dd));
305  if (dd > d) std::fill(this->begin() + alpha(n,d), this->end(), T(0));
306  d = dd;
307  }
308 
309  template<typename T>
310  void polynomial<T>::add_monomial(const T &coeff, const power_index &power) {
311  size_type i = power.global_index();
312  GMM_ASSERT2(n == power.size(), "dimensions mismatch");
313  if (i >= this->size()) { change_degree(power.degree()); }
314  ((*this)[i]) += coeff;
315  }
316 
317  template<typename T>
318  polynomial<T> &polynomial<T>::operator +=(const polynomial &Q) {
319  GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
320 
321  if (Q.degree() > degree()) change_degree(Q.degree());
322  iterator it = this->begin();
323  const_iterator itq = Q.begin(), ite = Q.end();
324  for ( ; itq != ite; ++itq, ++it) *it += *itq;
325  return *this;
326  }
327 
328  template<typename T>
329  polynomial<T> &polynomial<T>::operator -=(const polynomial &Q) {
330  GMM_ASSERT2(Q.dim() == dim() && dim() != 0, "dimensions mismatch");
331 
332  if (Q.degree() > degree()) change_degree(Q.degree());
333  iterator it = this->begin();
334  const_iterator itq = Q.begin(), ite = Q.end();
335  for ( ; itq != ite; ++itq, ++it) *it -= *itq;
336  return *this;
337  }
338 
339  template<typename T>
340  polynomial<T> polynomial<T>::operator -() const {
341  polynomial<T> Q = *this;
342  iterator itq = Q.begin(), ite = Q.end();
343  for ( ; itq != ite; ++itq) *itq = -(*itq);
344  return Q;
345  }
346 
347  template<typename T>
348  polynomial<T> &polynomial<T>::operator *=(const polynomial &Q) {
349  GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
350 
351  polynomial aux = *this;
352  change_degree(0); (*this)[0] = T(0);
353 
354  power_index miq(Q.dim()), mia(dim()), mitot(dim());
355  if (dim() > 0) miq[dim()-1] = Q.degree();
356  const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
357  for ( ; itq != ite; ++itq, --miq) {
358  if (*itq != T(0)) {
359  reverse_iterator ita = aux.rbegin(), itae = aux.rend();
360  std::fill(mia.begin(), mia.end(), 0);
361  if (dim() > 0) mia[dim()-1] = aux.degree();
362  for ( ; ita != itae; ++ita, --mia)
363  if (*ita != T(0)) {
364  power_index::iterator mita = mia.begin(), mitq = miq.begin();
365  power_index::iterator mit = mitot.begin(), mite = mia.end();
366  for ( ; mita != mite; ++mita, ++mitq, ++mit)
367  *mit = short_type((*mita) + (*mitq)); /* on pourrait calculer
368  directement l'index global. */
369  // cerr << "*= : " << *this << ", itq*ita="
370  // << (*itq) * (*ita) << endl;
371  // cerr << " itq = " << *itq << endl;
372  // cerr << " ita = " << *ita << endl;
373  add_monomial((*itq) * (*ita), mitot);
374 
375  }
376  }
377  }
378  return *this;
379  }
380 
381  template<typename T>
382  void polynomial<T>::direct_product(const polynomial &Q) {
383  polynomial aux = *this;
384 
385  change_degree(0); n = short_type(n+Q.dim()); (*this)[0] = T(0);
386 
387  power_index miq(Q.dim()), mia(aux.dim()), mitot(dim());
388  if (Q.dim() > 0) miq[Q.dim()-1] = Q.degree();
389  const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
390  for ( ; itq != ite; ++itq, --miq)
391  if (*itq != T(0)) {
392  reverse_iterator ita = aux.rbegin(), itae = aux.rend();
393  std::fill(mia.begin(), mia.end(), 0);
394  if (aux.dim() > 0) mia[aux.dim()-1] = aux.degree();
395  for ( ; ita != itae; ++ita, --mia)
396  if (*ita != T(0)) {
397  std::copy(mia.begin(), mia.end(), mitot.begin());
398  std::copy(miq.begin(), miq.end(), mitot.begin() + aux.dim());
399  add_monomial((*itq) * (*ita), mitot); /* on pourrait calculer
400  directement l'index global. */
401  }
402  }
403  }
404 
405  template<typename T>
406  polynomial<T> &polynomial<T>::operator *=(const T &e) {
407  iterator it = this->begin(), ite = this->end();
408  for ( ; it != ite; ++it) (*it) *= e;
409  return *this;
410  }
411 
412  template<typename T>
413  polynomial<T> &polynomial<T>::operator /=(const T &e) {
414  iterator it = this->begin(), ite = this->end();
415  for ( ; it != ite; ++it) (*it) /= e;
416  return *this;
417  }
418 
419  template<typename T>
420  inline void polynomial<T>::derivative(short_type k) {
421  GMM_ASSERT2(k < n, "index out of range");
422 
423  iterator it = this->begin(), ite = this->end();
424  power_index mi(dim());
425  for ( ; it != ite; ++it) {
426  if ((*it) != T(0) && mi[k] > 0)
427  { mi[k]--; (*this)[mi.global_index()] = (*it) * T(mi[k] + 1); mi[k]++; }
428  *it = T(0);
429  ++mi;
430  }
431  if (d > 0) change_degree(short_type(d-1));
432  }
433 
434  template<typename T> template<typename ITER>
435  inline T polynomial<T>::horner(power_index &mi, short_type k,
436  short_type de, const ITER &it) const {
437  if (k == 0)
438  return (*this)[mi.global_index()];
439  else {
440  T v = (*(it+k-1)), res = T(0);
441  for (mi[k-1] = short_type(degree() - de); mi[k-1] != short_type(-1);
442  (mi[k-1])--)
443  res = horner(mi, short_type(k-1), short_type(de + mi[k-1]), it)
444  + v * res;
445  mi[k-1] = 0;
446  return res;
447  }
448  }
449 
450 
451  template<typename T> template<typename ITER>
452  T polynomial<T>::eval(const ITER &it) const {
453  /* direct evaluation for common low degree polynomials */
454  unsigned deg = degree();
455  const_iterator P = this->begin();
456  if (deg == 0) return P[0];
457  else if (deg == 1) {
458  T s = P[0];
459  for (size_type i=0; i < dim(); ++i) s += it[i]*P[i+1];
460  return s;
461  }
462 
463  switch (dim()) {
464  case 1: {
465  T x = it[0];
466  if (deg == 2) return P[0] + x*(P[1] + x*(P[2]));
467  if (deg == 3) return P[0] + x*(P[1] + x*(P[2] + x*(P[3])));
468  if (deg == 4) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4]))));
469  if (deg == 5) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] + x*(P[5])))));
470  if (deg == 6) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] + x*(P[5] + x*(P[6]))))));
471  } break;
472  case 2: {
473  T x = it[0];
474  T y = it[1];
475  if (deg == 2) return P[0] + x*(P[1] + x*(P[3])) + y*(P[2] + x*(P[4]) + y*(P[5]));
476  if (deg == 3) return P[0] + x*(P[1] + x*(P[3] + x*(P[6]))) + y*(P[2] + x*(P[4] + x*(P[7])) + y*(P[5] + x*(P[8]) + y*(P[9])));
477  if (deg == 4) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10])))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11]))) + y*(P[5] + x*(P[8] + x*(P[12])) + y*(P[9] + x*(P[13]) + y*(P[14]))));
478  if (deg == 5) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] + x*(P[15]))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16])))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17]))) + y*(P[9] + x*(P[13] + x*(P[18])) + y*(P[14] + x*(P[19]) + y*(P[20])))));
479  if (deg == 6) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10] + x*(P[15] + x*(P[21])))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16] + x*(P[22]))))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17] + x*(P[23])))) + y*(P[9] + x*(P[13] + x*(P[18] + x*(P[24]))) + y*(P[14] + x*(P[19] + x*(P[25])) + y*(P[20] + x*(P[26]) + y*(P[27]))))));
480  } break;
481  case 3: {
482  T x = it[0];
483  T y = it[1];
484  T z = it[2];
485  if (deg == 2) return P[0] + x*(P[1] + x*(P[4])) + y*(P[2] + x*(P[5]) + y*(P[7])) + z*(P[3] + x*(P[6]) + y*(P[8]) + z*(P[9]));
486  if (deg == 3) return P[0] + x*(P[1] + x*(P[4] + x*(P[10]))) + y*(P[2] + x*(P[5] + x*(P[11])) + y*(P[7] + x*(P[13]) + y*(P[16]))) + z*(P[3] + x*(P[6] + x*(P[12])) + y*(P[8] + x*(P[14]) + y*(P[17])) + z*(P[9] + x*(P[15]) + y*(P[18]) + z*(P[19])));
487  if (deg == 4) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20])))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]))) + y*(P[7] + x*(P[13] + x*(P[23])) + y*(P[16] + x*(P[26]) + y*(P[30])))) + z*(P[3] + x*(P[6] + x*(P[12] + x*(P[22]))) + y*(P[8] + x*(P[14] + x*(P[24])) + y*(P[17] + x*(P[27]) + y*(P[31]))) + z*(P[9] + x*(P[15] + x*(P[25])) + y*(P[18] + x*(P[28]) + y*(P[32])) + z*(P[19] + x*(P[29]) + y*(P[33]) + z*(P[34]))));
488  if (deg == 5) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20] + x*(P[35]))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + x*(P[36])))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38]))) + y*(P[16] + x*(P[26] + x*(P[41])) + y*(P[30] + x*(P[45]) + y*(P[50]))))) + z*(P[3] + x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37])))) + y*(P[8] + x*(P[14] + x*(P[24] + x*(P[39]))) + y*(P[17] + x*(P[27] + x*(P[42])) + y*(P[31] + x*(P[46]) + y*(P[51])))) + z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40]))) + y*(P[18] + x*(P[28] + x*(P[43])) + y*(P[32] + x*(P[47]) + y*(P[52]))) + z*(P[19] + x*(P[29] + x*(P[44])) + y*(P[33] + x*(P[48]) + y*(P[53])) + z*(P[34] + x*(P[49]) + y*(P[54]) + z*(P[55])))));
489  if (deg == 6) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20] + x*(P[35] + x*(P[56])))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + x*(P[36] + x*(P[57]))))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38] + x*(P[59])))) + y*(P[16] + x*(P[26] + x*(P[41] + x*(P[62]))) + y*(P[30] + x*(P[45] + x*(P[66])) + y*(P[50] + x*(P[71]) + y*(P[77])))))) + z*(P[3] + x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37] + x*(P[58]))))) + y*(P[8] + x*(P[14] + x*(P[24] + x*(P[39] + x*(P[60])))) + y*(P[17] + x*(P[27] + x*(P[42] + x*(P[63]))) + y*(P[31] + x*(P[46] + x*(P[67])) + y*(P[51] + x*(P[72]) + y*(P[78]))))) + z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40] + x*(P[61])))) + y*(P[18] + x*(P[28] + x*(P[43] + x*(P[64]))) + y*(P[32] + x*(P[47] + x*(P[68])) + y*(P[52] + x*(P[73]) + y*(P[79])))) + z*(P[19] + x*(P[29] + x*(P[44] + x*(P[65]))) + y*(P[33] + x*(P[48] + x*(P[69])) + y*(P[53] + x*(P[74]) + y*(P[80]))) + z*(P[34] + x*(P[49] + x*(P[70])) + y*(P[54] + x*(P[75]) + y*(P[81])) + z*(P[55] + x*(P[76]) + y*(P[82]) + z*(P[83]))))));
490  } break;
491  }
492 
493  /*
494  switch (deg) {
495  case 0: return (*this)[0];
496  case 1: {
497  T s = (*this)[0];
498  for (size_type i=0; i < dim(); ++i) s += it[i]*(*this)[i+1];
499  return s;
500  }
501  case 2:
502  case 3: {
503  if (dim() == 1) {
504  const T &x = it[0];
505  if (deg == 2) return p[0] + x*(p[1] + x*p[2]);
506  else if (deg == 3) return p[0] + x*(p[1] + x*(p[2]+x*p[3]));
507  } else if (dim() == 2) {
508  const T &x = it[0];
509  const T &y = it[1];
510  if (deg == 2)
511  return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y;
512  else if (deg == 3)
513  return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y +
514  p[6]*x*x*x + p[7]*x*x*y + p[8]*x*y*y + p[9]*y*y*y;
515  } else if (dim() == 3) {
516  const T &x = it[0];
517  const T &y = it[1];
518  const T &z = it[2];
519  if (deg == 2)
520  return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + p[6]*x*z +
521  p[7]*y*y + p[8]*y*z + p[9]*z*z;
522  else if (deg == 3)
523  return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y + p[6]*x*z +
524  p[7]*y*y + p[8]*y*z + p[9]*z*z +
525  p[10]*x*x*x + p[11]*x*x*y + p[12]*x*x*z + p[13]*x*y*y + p[14]*x*y*z + p[15]*x*z*z +
526  p[16]*y*y*y + p[17]*y*y*z + p[18]*y*z*z +
527  p[19]*z*z*z;
528  }
529  }
530  }*/
531  /* for other polynomials, Horner evaluation (quite slow..) */
532  power_index mi(dim());
533  return horner(mi, dim(), 0, it);
534  }
535 
536  template<typename ITER>
537  typename std::iterator_traits<ITER>::value_type
538  eval_monomial(const power_index &mi, ITER it) {
539  typename std::iterator_traits<ITER>::value_type res
540  = typename std::iterator_traits<ITER>::value_type(1);
541  power_index::const_iterator mit = mi.begin(), mite = mi.end();
542  for ( ; mit != mite; ++mit, ++it)
543  for (short_type l = 0; l < *mit; ++l)
544  res *= *it;
545  return res;
546  }
547 
548 
549  /// Print P to the output stream o. for instance cout << P;
550  template<typename T> std::ostream &operator <<(std::ostream &o,
551  const polynomial<T>& P) {
552  bool first = true; size_type n = 0;
553  typename polynomial<T>::const_iterator it = P.begin(), ite = P.end();
554  power_index mi(P.dim());
555  if (it != ite && *it != T(0))
556  { o << *it; first = false; ++it; ++n; ++mi; }
557  for ( ; it != ite ; ++it, ++mi ) {
558  if (*it != T(0)) {
559  bool first_var = true;
560  if (!first) { if (*it < T(0)) o << " - "; else o << " + "; }
561  else if (*it < T(0)) o << "-";
562  if (gmm::abs(*it)!=T(1)) { o << gmm::abs(*it); first_var = false; }
563  for (short_type j = 0; j < P.dim(); ++j)
564  if (mi[j] != 0) {
565  if (!first_var) o << "*";
566  first_var = false;
567  if (P.dim() <= 7) o << "xyzwvut"[j];
568  else o << "x_" << j;
569  if (mi[j] > 1) o << "^" << mi[j];
570  }
571  first = false; ++n;
572  }
573  }
574  if (n == 0) o << "0";
575  return o;
576  }
577 
578  /**
579  polynomial variable substitution
580  @param P the original polynomial
581  @param S the substitution poly (not a multivariate one)
582  @param subs_dim : which variable is substituted
583  example: poly_subs(x+y*x^2, x+1, 0) = x+1 + y*(x+1)^2
584  */
585  template<typename T>
586  polynomial<T> poly_substitute_var(const polynomial<T>& P,
587  const polynomial<T>& S,
588  size_type subs_dim) {
589  GMM_ASSERT2(S.dim() == 1 && subs_dim < P.dim(),
590  "wrong arguments for polynomial substitution");
591  polynomial<T> res(P.dim(),0);
592  bgeot::power_index pi(P.dim());
593  std::vector< polynomial<T> > Spow(1);
594  Spow[0] = polynomial<T>(1, 0); Spow[0].one(); // Spow stores powers of S
595  for (size_type k=0; k < P.size(); ++k, ++pi) {
596  if (P[k] == T(0)) continue;
597  while (pi[subs_dim] >= Spow.size())
598  Spow.push_back(S*Spow.back());
599  const polynomial<T>& p = Spow[pi[subs_dim]];
600  bgeot::power_index pi2(pi);
601  for (short_type i=0; i < p.size(); ++i) {
602  pi2[subs_dim] = i;
603  res.add_monomial(p[i]*P[k],pi2);
604  }
605  }
606  return res;
607  }
608 
609  template<typename U, typename T>
610  polynomial<T> operator *(T c, const polynomial<T> &p)
611  { polynomial<T> q = p; q *= c; return q; }
612 
613  typedef polynomial<opt_long_scalar_type> base_poly;
614 
615  /* usual constant polynomials */
616 
617  inline base_poly null_poly(short_type n) { return base_poly(n, 0); }
618  inline base_poly one_poly(short_type n)
619  { base_poly res=base_poly(n, 0); res.one(); return res; }
620  inline base_poly one_var_poly(short_type n, short_type k)
621  { return base_poly(n, 1, k); }
622 
623  /** read a base_poly on the stream ist. */
624  base_poly read_base_poly(short_type n, std::istream &f);
625 
626  /** read a base_poly on the string s. */
627  base_poly read_base_poly(short_type n, const std::string &s);
628 
629 
630  /**********************************************************************/
631  /* A class for rational fractions */
632  /**********************************************************************/
633 
634  template<typename T> class rational_fraction : public std::vector<T> {
635  protected :
636 
637  polynomial<T> numerator_, denominator_;
638 
639  public :
640 
641  const polynomial<T> &numerator() const { return numerator_; }
642  const polynomial<T> &denominator() const { return denominator_; }
643 
644  short_type dim() const { return numerator_.dim(); }
645 
646  /// Add Q to P. P contains the result.
647  rational_fraction &operator +=(const rational_fraction &Q) {
648  numerator_ = numerator_*Q.denominator() + Q.numerator()*denominator_;
649  denominator_ *= Q.denominator();
650  return *this;
651  }
652  /// Subtract Q from P. P contains the result.
653  rational_fraction &operator -=(const rational_fraction &Q) {
654  numerator_ = numerator_*Q.denominator() - Q.numerator()*denominator_;
655  denominator_ *= Q.denominator();
656  return *this;
657  }
658  /// Add Q to P.
659  rational_fraction operator +(const rational_fraction &Q) const
660  { rational_fraction R = *this; R += Q; return R; }
661  /// Subtract Q from P.
662  rational_fraction operator -(const rational_fraction &Q) const
663  { rational_fraction R = *this; R -= Q; return R; }
664  /// Add Q to P.
665  rational_fraction operator +(const polynomial<T> &Q) const
666  { rational_fraction R(numerator_+Q*denominator_, denominator_); return R; }
667  /// Subtract Q from P.
668  rational_fraction operator -(const polynomial<T> &Q) const
669  { rational_fraction R(numerator_-Q*denominator_, denominator_); return R; }
670  rational_fraction operator -() const
671  { rational_fraction R(-numerator_, denominator_); return R; }
672  /// Multiply P with Q. P contains the result.
673  rational_fraction &operator *=(const rational_fraction &Q)
674  { numerator_*=Q.numerator(); denominator_*=Q.denominator(); return *this; }
675  /// Divide P by Q. P contains the result.
676  rational_fraction &operator /=(const rational_fraction &Q)
677  { numerator_*=Q.denominator(); denominator_*=Q.numerator(); return *this; }
678  /// Multiply P with Q.
679  rational_fraction operator *(const rational_fraction &Q) const
680  { rational_fraction R = *this; R *= Q; return R; }
681  /// Divide P by Q.
682  rational_fraction operator /(const rational_fraction &Q) const
683  { rational_fraction R = *this; R /= Q; return R; }
684  /// Multiply P with the scalar a. P contains the result.
685  rational_fraction &operator *=(const T &e)
686  { numerator_ *= e; return *this; }
687  /// Multiply P with the scalar a.
688  rational_fraction operator *(const T &e) const
689  { rational_fraction R = *this; R *= e; return R; }
690  /// Divide P with the scalar a. P contains the result.
691  rational_fraction &operator /=(const T &e)
692  { denominator_ *= e; return *this; }
693  /// Divide P with the scalar a.
694  rational_fraction operator /(const T &e) const
695  { rational_fraction res = *this; res /= e; return res; }
696  /// operator ==.
697  bool operator ==(const rational_fraction &Q) const
698  { return (numerator_==Q.numerator() && denominator_==Q.denominator()); }
699  /// operator !=.
700  bool operator !=(const rational_fraction &Q) const
701  { return !(operator ==(*this,Q)); }
702  /// Derivative of P with respect to the variable k. P contains the result.
703  void derivative(short_type k) {
704  polynomial<T> der_num = numerator_; der_num.derivative(k);
705  polynomial<T> der_den = denominator_; der_den.derivative(k);
706  if (der_den.is_zero()) {
707  if (der_num.is_zero()) this->clear();
708  else numerator_ = der_num;
709  } else {
710  numerator_ = der_num * denominator_ - der_den * numerator_;
711  denominator_ = denominator_ * denominator_;
712  }
713  }
714  /// Makes P = 1.
715  void one() { numerator_.one(); denominator_.one(); }
716  void clear() { numerator_.clear(); denominator_.one(); }
717  template <typename ITER> T eval(const ITER &it) const {
718  typedef typename gmm::number_traits<T>::magnitude_type R;
719  T a = numerator_.eval(it), b = denominator_.eval(it);
720  if (b == T(0)) { // The better should be to evaluate the derivatives ...
721  std::vector<T> p(it, it+dim());
722  R no = gmm::vect_norm2(p);
723  if (no == R(0)) { gmm::fill_random(p); gmm::scale(p, R(1E-35)); }
724  else gmm::scale(p, R(0.9999999));
725  a = numerator_.eval(p.begin());
726  b = denominator_.eval(p.begin());
727  }
728  if (a != T(0)) a /= b;
729  return a;
730  }
731  /// Constructor.
732  rational_fraction()
733  : numerator_(1,0), denominator_(1,0) { clear(); }
734  /// Constructor.
735  rational_fraction(short_type dim_)
736  : numerator_(dim_,0), denominator_(dim_,0) { clear(); }
737  /// Constructor
738  rational_fraction(const polynomial<T> &numer)
739  : numerator_(numer), denominator_(numer.dim(),0) { denominator_.one(); }
740  /// Constructor
741  rational_fraction(const polynomial<T> &numer, const polynomial<T> &denom)
742  : numerator_(numer), denominator_(denom)
743  { GMM_ASSERT1(numer.dim() == denom.dim(), "Dimensions mismatch"); }
744  };
745 
746  /// Add Q to P.
747  template<typename T>
748  rational_fraction<T> operator +(const polynomial<T> &P,
749  const rational_fraction<T> &Q) {
750  rational_fraction<T> R(P*Q.denominator()+Q.numerator(), Q.denominator());
751  return R;
752  }
753  /// Subtract Q from P.
754  template<typename T>
755  rational_fraction<T> operator -(const polynomial<T> &P,
756  const rational_fraction<T> &Q) {
757  rational_fraction<T> R(P*Q.denominator()-Q.numerator(), Q.denominator());
758  return R;
759  }
760 
761  template<typename T> std::ostream &operator <<
762  (std::ostream &o, const rational_fraction<T>& P) {
763  o << "[" << P.numerator() << "]/[" << P.denominator() << "]";
764  return o;
765  }
766 
767  typedef rational_fraction<opt_long_scalar_type> base_rational_fraction;
768 
769 } /* end of namespace bgeot. */
770 
771 
772 #endif /* BGEOT_POLY_H__ */
bgeot_config.h
defines and typedefs for namespace bgeot
bgeot::polynomial
This class deals with plain polynomials with several variables.
Definition: bgeot_poly.h:179
bgeot::polynomial::real_degree
short_type real_degree() const
gives the degree of the polynomial, considering only non-zero coefficients
Definition: bgeot_poly.h:294
bgeot::polynomial::polynomial
polynomial(short_type dim_, short_type degree_, short_type k)
Constructor for the polynomial 'x' (k=0), 'y' (k=1), 'z' (k=2) etc.
Definition: bgeot_poly.h:266
bgeot::polynomial::operator+=
polynomial & operator+=(const polynomial &Q)
Add Q to P. P contains the result.
Definition: bgeot_poly.h:318
bgeot::polynomial::dim
short_type dim() const
Gives the dimension (number of variables)
Definition: bgeot_poly.h:198
bgeot::polynomial::operator-=
polynomial & operator-=(const polynomial &Q)
Subtract Q from P. P contains the result.
Definition: bgeot_poly.h:329
bgeot::polynomial::polynomial
polynomial()
Constructor.
Definition: bgeot_poly.h:253
bgeot::polynomial::add_monomial
void add_monomial(const T &coeff, const power_index &power)
Add to the polynomial a monomial of coefficient a and correpsonding to the power index pi.
Definition: bgeot_poly.h:310
bgeot::polynomial::operator*
polynomial operator*(const polynomial &Q) const
Multiply P with Q.
Definition: bgeot_poly.h:275
bgeot::polynomial::polynomial
polynomial(short_type dim_, short_type degree_)
Constructor.
Definition: bgeot_poly.h:262
bgeot::polynomial::eval
T eval(const ITER &it) const
Evaluate the polynomial.
Definition: bgeot_poly.h:452
bgeot::polynomial::change_degree
void change_degree(short_type dd)
Change the degree of the polynomial to d.
Definition: bgeot_poly.h:303
bgeot::polynomial::derivative
void derivative(short_type k)
Derivative of P with respect to the variable k. P contains the result.
Definition: bgeot_poly.h:420
bgeot::polynomial::operator+
polynomial operator+(const polynomial &Q) const
Add Q to P.
Definition: bgeot_poly.h:210
bgeot::polynomial::operator!=
bool operator!=(const polynomial &Q) const
operator !=.
Definition: bgeot_poly.h:236
bgeot::polynomial::one
void one()
Makes P = 1.
Definition: bgeot_poly.h:241
bgeot::polynomial::operator*=
polynomial & operator*=(const polynomial &Q)
Multiply P with Q. P contains the result.
Definition: bgeot_poly.h:348
bgeot::polynomial::degree
short_type degree() const
Gives the degree of the polynomial.
Definition: bgeot_poly.h:192
bgeot::polynomial::operator/=
polynomial & operator/=(const T &e)
Divide P with the scalar a. P contains the result.
Definition: bgeot_poly.h:413
bgeot::polynomial::operator/
polynomial operator/(const T &e) const
Divide P with the scalar a.
Definition: bgeot_poly.h:231
bgeot::polynomial::direct_product
void direct_product(const polynomial &Q)
Product of P and Q considering that variables of Q come after variables of P.
Definition: bgeot_poly.h:382
bgeot::polynomial::operator==
bool operator==(const polynomial &Q) const
operator ==.
Definition: bgeot_poly.h:279
bgeot::power_index
Vector of integer (16 bits type) which represent the powers of a monomial.
Definition: bgeot_poly.h:63
bgeot::power_index::power_index
power_index()
Constructor.
Definition: bgeot_poly.h:107
bgeot::power_index::operator++
const power_index & operator++()
Gives the next power index.
Definition: bgeot_poly.cc:53
bgeot::power_index::degree
short_type degree() const
Gives the degree.
Definition: bgeot_poly.cc:89
bgeot::power_index::global_index
size_type global_index() const
Gives the global number of the index (i.e.
Definition: bgeot_poly.cc:95
bgeot::power_index::operator--
const power_index & operator--()
Gives the next previous index.
Definition: bgeot_poly.cc:70
dal_static_stored_objects.h
Stores interdependent getfem objects.
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:556
gmm::fill_random
void fill_random(L &l)
fill a vector or matrix with random value (uniform [-1,1]).
Definition: gmm_blas.h:129
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:26
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::read_base_poly
base_poly read_base_poly(short_type n, std::istream &f)
read a base_poly on the stream ist.
Definition: bgeot_poly.cc:210
bgeot::operator-
rational_fraction< T > operator-(const polynomial< T > &P, const rational_fraction< T > &Q)
Subtract Q from P.
Definition: bgeot_poly.h:755
bgeot::poly_substitute_var
polynomial< T > poly_substitute_var(const polynomial< T > &P, const polynomial< T > &S, size_type subs_dim)
polynomial variable substitution
Definition: bgeot_poly.h:586
bgeot::operator+
rational_fraction< T > operator+(const polynomial< T > &P, const rational_fraction< T > &Q)
Add Q to P.
Definition: bgeot_poly.h:748
bgeot::operator<<
std::ostream & operator<<(std::ostream &o, const convex_structure &cv)
Print the details of the convex structure cvs to the output stream o.
Definition: bgeot_convex_structure.cc:70
bgeot::operator==
bool operator==(const pconvex_structure &p1, const pconvex_structure &p2)
Stored objects must be compared by keys, because there is a possibility that they are duplicated in s...
Definition: bgeot_convex_structure.cc:110
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::alpha
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
Definition: bgeot_poly.cc:46