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