GetFEM  5.4.2
bgeot_poly.cc
1 /*===========================================================================
2 
3  Copyright (C) 2000-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 
23 
24 #include "getfem/bgeot_poly.h"
25 #include "getfem/bgeot_small_vector.h"
26 #include "getfem/bgeot_ftool.h"
27 #include "getfem/getfem_locale.h"
28 
29 namespace bgeot {
30 
31 #define STORED 150
32  static gmm::dense_matrix<size_type> alpha_M_(STORED, STORED);
33  static void alpha_init_() {
34  static bool init = false;
35  if (!init) {
36  for (short_type i = 0; i < STORED; ++i) {
37  alpha_M_(i, 0) = alpha_M_(0, i) = 1;
38  for (short_type j = 1; j <= i; ++j)
39  alpha_M_(i,j) = alpha_M_(j,i) = (alpha_M_(i, j-1) * (i+j)) / j;
40  }
41  init = true;
42  }
43  }
44  static inline size_type alpha_(short_type n, short_type d)
45  { return alpha_M_(d,n); }
46 
47  size_type alpha(short_type n, short_type d) {
48  alpha_init_();
49  GMM_ASSERT1(n < STORED && d < STORED,
50  "alpha called with n = " << n << " and d = " << d);
51  return alpha_(n,d);
52  }
53 
54  const power_index &power_index::operator ++() {
55  short_type n = short_type(size()), l;
56  if (n > 0) {
57  size_type g_idx = global_index_; short_type deg = degree_;
58  reverse_iterator it = rbegin() + 1;
59  for (l = short_type(n-2); l != short_type(-1); --l, ++it)
60  if (*it != 0) break;
61  short_type a = (*this)[n-1]; (*this)[n-1] = 0;
62  (*this)[short_type(l+1)] = short_type(a + 1);
63  if (l != short_type(-1)) ((*this)[l])--;
64  else if (short_type(deg+1)) degree_ = short_type(deg+1);
65  if (g_idx+1) global_index_ = g_idx+1;
66  //degree_ = short_type(-1);
67  }
68  return *this;
69  }
70 
71  const power_index &power_index::operator --() {
72  short_type n = short_type(size()), l;
73  if (n > 0) {
74  size_type g_idx = global_index_; short_type deg = degree_;
75  reverse_iterator it = rbegin();
76  for (l = short_type(n-1); l != short_type(-1); --l, ++it) {
77  if (*it != 0) break;
78  }
79  if (l != short_type(-1)) {
80  short_type a = (*this)[l];
81  (*this)[l] = 0; (*this)[n-1] = short_type(a - 1);
82  if (l > 0) ((*this)[l-1])++;
83  else if (short_type(deg+1)) degree_ = short_type(deg-1);
84  }
85  if (g_idx+1) global_index_ = g_idx-1;
86  }
87  return *this;
88  }
89 
90  short_type power_index::degree() const {
91  if (degree_ != short_type(-1)) return degree_;
92  degree_ = short_type(std::accumulate(begin(), end(), 0));
93  return degree_;
94  }
95 
96  size_type power_index::global_index(void) const {
97  if (global_index_ != size_type(-1)) return global_index_;
98  short_type d = degree(), n = short_type(size());
99  global_index_ = 0;
100  const_iterator it = begin(), ite = end();
101  for ( ; it != ite && d > 0; ++it)
102  { global_index_ += alpha_(n,short_type(d-1)); d=short_type(d-*it); --n; }
103  return global_index_;
104  }
105 
106  power_index::power_index(short_type nn) : v(nn), degree_(0), global_index_(0)
107  { std::fill(begin(), end(), short_type(0)); alpha_init_(); }
108 
109 
110  // functions to read a polynomial on a stream
111 
112  static void parse_error(int i)
113  { GMM_ASSERT1(false, "Syntax error reading a polynomial " << i); }
114 
115  static std::string stored_s;
116  int stored_tokent;
117 
118  static void unget_token(int i, std::string s)
119  { std::swap(s, stored_s); stored_tokent = i; }
120 
121  static int get_next_token(std::string &s, std::istream &f) {
122  if (stored_s.size() == 0) {
123  int r = get_token(f, s, true, false, false);
124  return r;
125  }
126  else { s.clear(); std::swap(s, stored_s); return stored_tokent; }
127  }
128 
129  static base_poly read_expression(short_type n, std::istream &f) {
130  gmm::stream_standard_locale sl(f);
131  getfem::standard_locale sll;
132  base_poly result(n,0);
133  std::string s;
134  int i = get_next_token(s, f), j;
135  switch (i) {
136  case 2 : result.one();
137  result *= opt_long_scalar_type(::strtod(s.c_str(), 0));
138  break;
139  case 4 :
140  if (s == "x") result = base_poly(n, 1, 0);
141  else if (s == "y" && n > 1) result = base_poly(n, 1, 1);
142  else if (s == "z" && n > 2) result = base_poly(n, 1, 2);
143  else if (s == "w" && n > 3) result = base_poly(n, 1, 3);
144  else if (s == "v" && n > 4) result = base_poly(n, 1, 4);
145  else if (s == "u" && n > 5) result = base_poly(n, 1, 5);
146  else if (s == "t" && n > 6) result = base_poly(n, 1, 6);
147  else if (s == "sqrt") {
148  base_poly p = read_expression(n, f);
149  if (p.degree() > 0) parse_error(1);
150  result.one(); result *= sqrt(p[0]);
151  }
152  else { parse_error(2); }
153  break;
154  case 5 :
155  switch (s[0]) {
156  case '(' :
157  result = read_base_poly(n, f);
158  j = get_next_token(s, f);
159  if (j != 5 || s[0] != ')') parse_error(3);
160  break;
161  default : parse_error(4);
162  }
163  break;
164  default : parse_error(5);
165  }
166  return result;
167  }
168 
169  static void operator_priority_(int i, char c, int &prior, int &op) {
170  if (i == 5)
171  switch (c) {
172  case '*' : prior = 2; op = 1; return;
173  case '/' : prior = 2; op = 2; return;
174  case '+' : prior = 3; op = 3; return;
175  case '-' : prior = 3; op = 4; return;
176  case '^' : prior = 1; op = 5; return;
177  }
178  prior = op = 0;
179  }
180 
181  void do_bin_op(std::vector<base_poly> &value_list,
182  std::vector<int> &op_list,
183  std::vector<int> &prior_list) {
184  base_poly &p2 = *(value_list.rbegin());
185  if (op_list.back() != 6) {
186  assert(value_list.size()>1);
187  base_poly &p1 = *(value_list.rbegin()+1);
188  switch (op_list.back()) {
189  case 1 : p1 *= p2; break;
190  case 2 : if (p2.degree() > 0) parse_error(6); p1 /= p2[0]; break;
191  case 3 : p1 += p2; break;
192  case 4 : p1 -= p2; break;
193  case 5 :
194  {
195  if (p2.degree() > 0) parse_error(7);
196  int pow = int(to_scalar(p2[0]));
197  if (p2[0] != opt_long_scalar_type(pow) || pow < 0) parse_error(8);
198  base_poly p = p1; p1.one();
199  for (int i = 0; i < pow; ++i) p1 *= p;
200  }
201  break;
202  default: assert(0);
203  }
204  value_list.pop_back();
205  } else {
206  p2 *= opt_long_scalar_type(-1);
207  }
208  op_list.pop_back(); prior_list.pop_back();
209  }
210 
211  base_poly read_base_poly(short_type n, std::istream &f) {
212  std::vector<base_poly> value_list;
213  std::string s;
214  std::vector<int> op_list, prior_list;
215 
216  int i = get_next_token(s, f), prior, op;
217  if (i == 5 && s[0] == '-')
218  { op_list.push_back(6); prior_list.push_back(2); }
219  else if (i == 5 && s[0] == '+') ;
220  else unget_token(i, s);
221 
222  value_list.push_back(read_expression(n, f));
223  i = get_next_token(s, f);
224  operator_priority_(i, i ? s[0] : '0', prior, op);
225  while (op) {
226  while (!prior_list.empty() && prior_list.back() <= prior)
227  do_bin_op(value_list, op_list, prior_list);
228 
229  value_list.push_back(read_expression(n, f));
230  op_list.push_back(op);
231  prior_list.push_back(prior);
232 
233  i = get_next_token(s, f);
234  operator_priority_(i, i ? s[0] : '0', prior, op);
235  }
236 
237  if (i == 5 && s[0] == ')') { f.putback(')'); }
238  else if (i != 0 && (i != 5 || s[0] != ';')) {
239  cout << "s = " << s << endl;
240  parse_error(9);
241  }
242 
243  while (!prior_list.empty()) do_bin_op(value_list, op_list, prior_list);
244 
245  return value_list[0];
246  }
247 
248  base_poly read_base_poly(short_type n, const std::string &s) {
249  std::stringstream f(s);
250  return read_base_poly(n, f);
251  }
252 
253 
254 } /* end of namespace bgeot. */
bgeot::get_token
int get_token(std::istream &ist, std::string &st, bool ignore_cr, bool to_up, bool read_un_pm, int *linenb)
Very simple lexical analysis of general interest for reading small languages with a "MATLAB like" syn...
Definition: bgeot_ftool.cc:50
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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::power_index::degree
short_type degree() const
Gives the degree.
Definition: bgeot_poly.cc:90
getfem::standard_locale
Identical to gmm::standard_locale, but does not change std::locale in multi-threaded sections of the ...
Definition: getfem_locale.h:50
bgeot::power_index::operator--
const power_index & operator--()
Gives the next previous index.
Definition: bgeot_poly.cc:71
bgeot_poly.h
Multivariate polynomials.
getfem_locale.h
thread safe standard locale with RAII semantics
bgeot::polynomial< opt_long_scalar_type >
bgeot::power_index
Vector of integer (16 bits type) which represent the powers of a monomial.
Definition: bgeot_poly.h:64
bgeot_ftool.h
"File Tools"
bgeot::power_index::global_index
size_type global_index() const
Gives the global number of the index (i.e.
Definition: bgeot_poly.cc:96
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_small_vector.h
Small (dim < 8) vectors.
bgeot
Basic Geometric Tools.
Definition: bgeot_convex_ref.cc:27
bgeot::power_index::power_index
power_index()
Constructor.
Definition: bgeot_poly.h:108