Code-Eli  0.3.6
polynomial.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_mutil_poly_polynomial_hpp
14 #define eli_mutil_poly_polynomial_hpp
15 
16 #include <cmath>
17 #include <iterator>
18 #include <limits>
19 
20 #include "eli/code_eli.hpp"
21 
23 
24 namespace eli
25 {
26  namespace mutil
27  {
28  namespace poly
29  {
30  template<typename data__>
31  class polynomial
32  {
33  public:
34  typedef data__ data_type;
35  typedef Eigen::Matrix<data_type, Eigen::Dynamic, 1> coefficient_type;
36  typedef typename coefficient_type::Index index_type;
37 
38  private:
39  coefficient_type a;
40 
41  public:
43  polynomial(const coefficient_type &c) : a(c) {}
44  polynomial(const polynomial<data__> &p) : a(p.a) {}
45  polynomial(const data_type &root) {set_roots(root);}
46  polynomial(const data_type &root1, const data_type &root2) {set_roots(root1, root2);}
47  polynomial(const data_type &root1, const data_type &root2, const data_type &root3) {set_roots(root1, root2, root3);}
48  polynomial(const data_type &root1, const data_type &root2, const data_type &root3, const data_type &root4){set_roots(root1, root2, root3, root4);}
49  template<typename itroot__>
50  polynomial(itroot__ its, itroot__ ite) {set_roots(its, ite);}
51 
52  polynomial<data__> & operator=(const data_type &d)
53  {
54  coefficient_type a_new(1);
55 
56  a_new(0)=d;
57  set_coefficients(a_new);
58 
59  return *this;
60  }
61 
63  {
64  if (this!=&p)
65  {
66  a=p.a;
67  }
68 
69  return *this;
70  }
71 
72  index_type degree() const
73  {
74  return size()-1;
75  }
76 
77  data_type coefficient(const index_type &i) const
78  {
79  return a(i);
80  }
81  void set_coefficients(const coefficient_type &ain)
82  {
83  a=ain;
84  }
85 
86  void get_coefficients(coefficient_type &aout) const
87  {
88  aout=a;
89  }
90 
91  void compress()
92  {
93  coefficient_type a_new;
94  index_type i,j;
95  for (i=this->degree(); i>0; --i)
96  {
97  if (a(i)!=static_cast<data_type>(0))
98  break;
99  }
100 
101  a_new.resize(i+1);
102  for (j=0; j<=i; ++j)
103  a_new(j)=a(j);
104 
105  set_coefficients(a_new);
106  }
107 
108  void adjust_zero(const data_type &small)
109  {
110  for (index_type i=0; i<this->size(); ++i)
111  {
112  if (std::abs(a(i))<small)
113  a(i)=static_cast<data_type>(0);
114  }
115  compress();
116  }
117 
118  void set_roots(const data_type &root)
119  {
120  coefficient_type a_new(2);
121 
122  a_new(1)=1;
123  a_new(0)=-root;
124  set_coefficients(a_new);
125  }
126 
127  void set_roots(const data_type &root1, const data_type &root2)
128  {
129  coefficient_type a_new(3);
130 
131  a_new(2)=1;
132  a_new(1)=-(root1+root2);
133  a_new(0)=root1*root2;
134  set_coefficients(a_new);
135  }
136 
137  void set_roots(const data_type &root1, const data_type &root2, const data_type &root3)
138  {
139  coefficient_type a_new(4);
140 
141  a_new(3)=1;
142  a_new(2)=-(root1+root2+root3);
143  a_new(1)=root1*root2+root1*root3+root2*root3;
144  a_new(0)=-root1*root2*root3;
145  set_coefficients(a_new);
146  }
147 
148  void set_roots(const data_type &root1, const data_type &root2, const data_type &root3, const data_type &root4)
149  {
150  coefficient_type a_new(5);
151 
152  a_new(4)=1;
153  a_new(3)=-(root1+root2+root3+root4);
154  a_new(2)=root1*root2+root1*root3+root1*root4+root2*root3+root2*root4+root3*root4;
155  a_new(1)=-(root1*root2*root3+root1*root2*root4+root1*root3*root4+root2*root3*root4);
156  a_new(0)=root1*root2*root3*root4;
157  set_coefficients(a_new);
158  }
159 
164  template<typename itroot__>
165  void set_roots(itroot__ its, itroot__ ite)
166  {
167  int deg=static_cast<int>(std::distance(its, ite)), n=deg+1;
168  coefficient_type a_new(n);
169  std::vector<data_type> root(deg);
170  std::vector<int> index(deg);
171  int i, j;
172 
173  // fill vector of roots
174  i=0;
175  for (itroot__ it=its; it!=ite; ++it, ++i)
176  {
177  index[i]=i;
178  root[i]=(*it);
179  }
180 
181  a_new[deg]=1;
182  for (j=1; j<n; ++j)
183  {
184  a_new[deg-j]=static_cast<data_type>(0);
185  do
186  {
187  data_type term(root[index[0]]);
188 
189  for (i=1; i<j; ++i)
190  term*=root[index[i]];
191  a_new[deg-j]+=term;
192  } while(eli::mutil::dm::next_combination(index.begin(), index.begin()+j, index.end()));
193  if (j%2==1)
194  a_new[deg-j]*=-1;
195  }
196  set_coefficients(a_new);
197  }
198 
199  data_type f(const data_type &t) const
200  {
201  data_type rtn(0);
202 
203  // check to make sure have valid curve
204  if (this->size()<=0)
205  {
206  assert(false);
207  return rtn;
208  }
209 
210  index_type i, n(this->degree());
211  for (i=n; i>0; --i)
212  {
213  rtn=t*(a(i)+rtn);
214  }
215  i=0;
216  rtn+=a(i);
217 
218  return rtn;
219  }
220 
221  data_type fp(const data_type &t) const
222  {
223  data_type rtn(0);
224 
225  // check to make sure have valid curve
226  if (this->size()<=0)
227  {
228  assert(false);
229  return rtn;
230  }
231 
232  index_type i, n;
233  n=static_cast<index_type>(this->degree());
234  for (i=n; i>1; --i)
235  {
236  rtn=t*(static_cast<data_type>(i)*a(i)+rtn);
237  }
238  i=1;
239  rtn+=static_cast<data_type>(i)*a(i);
240 
241  return rtn;
242  }
243 
244  data_type fpp(const data_type &t) const
245  {
246  data_type rtn(0);
247 
248  // check to make sure have valid curve
249  if (this->size()<=0)
250  {
251  assert(false);
252  return rtn;
253  }
254 
255  index_type i, n;
256  n=static_cast<index_type>(this->degree());
257  for (i=n; i>2; --i)
258  {
259  rtn=t*(static_cast<data_type>(i)*static_cast<data_type>(i-1)*a(i)+rtn);
260  }
261  i=2;
262  rtn+=static_cast<data_type>(i)*static_cast<data_type>(i-1)*a(i);
263 
264  return rtn;
265  }
266 
267  data_type fppp(const data_type &t) const
268  {
269  data_type rtn(0);
270 
271  // check to make sure have valid curve
272  if (this->size()<=0)
273  {
274  assert(false);
275  return rtn;
276  }
277 
278  index_type i, n;
279  n=static_cast<index_type>(this->degree());
280  for (i=n; i>3; --i)
281  {
282  rtn=t*(static_cast<data_type>(i)*static_cast<data_type>(i-1)*static_cast<data_type>(i-2)*a(i)+rtn);
283  }
284  i=3;
285  rtn+=static_cast<data_type>(i)*static_cast<data_type>(i-1)*static_cast<data_type>(i-2)*a(i);
286 
287  return rtn;
288  }
289 
291  {
292  if (this->size()<=0)
293  return 0;
294 
295  return new polynomial<data_type>(a);
296  }
297 
299  {
300  if (this->size()<=0)
301  return 0;
302 
303  index_type new_deg(this->degree()-1), deg(this->degree());
304  coefficient_type a_new(new_deg+1);
305  for (index_type i=1; i<=deg; ++i)
306  a_new(i-1)=static_cast<data_type>(i)*a(i);
307 
308  return new polynomial<data_type>(a_new);
309  }
310 
312  {
313  if (this->size()<=0)
314  return 0;
315 
316  index_type new_deg(this->degree()-2), deg(this->degree());
317  coefficient_type a_new(new_deg+1);
318  for (index_type i=2; i<=deg; ++i)
319  a_new(i-2)=static_cast<data_type>(i)*static_cast<data_type>(i-1)*a(i);
320 
321  return new polynomial<data_type>(a_new);
322  }
323 
325  {
326  if (this->size()<=0)
327  return 0;
328 
329  index_type new_deg(this->degree()-3), deg(this->degree());
330  coefficient_type a_new(new_deg+1);
331  for (index_type i=3; i<=deg; ++i)
332  a_new(i-3)=static_cast<data_type>(i)*static_cast<data_type>(i-1)*static_cast<data_type>(i-2)*a(i);
333 
334  return new polynomial<data_type>(a_new);
335  }
336 
338  {
339  // resize coefficients
340  index_type deg1(p1.degree()), deg2(p2.degree()), deg, n, i;
341 
342  deg=std::max(deg1, deg2);
343  n=deg+1;
344  coefficient_type a_new(n);
345  for (i=0; i<n; ++i)
346  {
347  if (i<=deg1)
348  a_new(i)=p1.a(i);
349  else
350  a_new(i)=static_cast<data_type>(0);
351 
352  if (i<=deg2)
353  a_new(i)+=p2.a(i);
354  }
355 
356  // set new coefficients
357  set_coefficients(a_new);
358  }
359 
361  {
362  // resize coefficients
363  index_type deg1(p1.degree()), deg2(p2.degree()), deg, n, i;
364 
365  deg=std::max(deg1, deg2);
366  n=deg+1;
367  coefficient_type a_new(n);
368  for (i=0; i<n; ++i)
369  {
370  if (i<=deg1)
371  a_new(i)=p1.a(i);
372  else
373  a_new(i)=static_cast<data_type>(0);
374 
375  if (i<=deg2)
376  a_new(i)-=p2.a(i);
377  }
378 
379  // set new coefficients
380  set_coefficients(a_new);
381  }
382 
384  {
385  // resize coefficients
386  index_type deg1(p1.degree()), deg2(p2.degree()), deg, n, j1, j2;
387 
388  deg=deg1+deg2;
389  n=deg+1;
390  coefficient_type a_new(n);
391  a_new.setZero();
392  for (j1=0; j1<=deg1; ++j1)
393  {
394  for (j2=0; j2<=deg2; ++j2)
395  {
396  a_new(j1+j2)+=p1.a(j1)*p2.a(j2);
397  }
398  }
399 
400  // set new coefficients
401  set_coefficients(a_new);
402  }
403 
404  void multiply(const data_type &d)
405  {
406  for (index_type i=0; i<this->size(); ++i)
407  a[i]*=d;
408  }
409 
411  {
412  // resize coefficients
413  index_type deg1(p1.degree()), deg2(p2.degree()), deg, n, j, j1, j2, nz;
414 
415  // catch case when deg2>deg1
416  if (deg2>deg1)
417  {
418  prem=p1;
419  coefficient_type a_new(1);
420  a_new.setZero();
421  set_coefficients(a_new);
422  return;
423  }
424 
425  // find the number of zero terms in the divisor to set the sizes
426  for (j=0; j<=deg2; ++j)
427  {
428  if (p2.a(deg2-j)!=static_cast<data_type>(0))
429  break;
430  }
431  nz=j;
432 
433  // if p2=0 then set results to nan and be done
434  if(nz>deg2)
435  {
436  coefficient_type a_new(1);
437  a_new(1)=std::numeric_limits<data__>::quiet_NaN();
438  prem.set_coefficients(a_new);
439  set_coefficients(a_new);
440  return;
441  }
442 
443  // set the size of quotient coefficients
444  deg=deg1+nz-deg2;
445  n=deg+1;
446  coefficient_type a_new(n);
447  a_new.setZero();
448 
449  // store current polynomial into remainder polynomial
450  polynomial<data_type> p1_copy(p1);
451 
452  // perform the division
453  deg2-=nz;
454 
455  for (j=0; j<=deg; ++j)
456  {
457  j1=deg-j;
458  a_new(j1)=p1_copy.a(deg2+j1)/p2.a(deg2);
459  for (j2=0; j2<=deg2; ++j2)
460  {
461  p1_copy.a(j2+j1)-=a_new(j1)*p2.a(j2);
462  }
463  }
464 
465  index_type nrem(deg2+nz);
466  coefficient_type a_rem(nrem);
467  for (j=0; j<nrem; ++j)
468  a_rem(j)=p1_copy.a(j);
469 
470  // set new coefficients
471  prem.set_coefficients(a_rem);
472  set_coefficients(a_new);
473  }
474 
475  void divide(const data_type &d)
476  {
477  for (index_type i=0; i<this->size(); ++i)
478  a[i]/=d;
479  }
480 
481  void negative()
482  {
483  for (index_type i=0; i<this->size(); ++i)
484  a[i]*=static_cast<data_type>(-1);
485  }
486 
487  private:
488  index_type size() const
489  {
490  return a.rows();
491  }
492 
493  };
494  }
495  }
496 }
497 
498 #endif
void set_roots(itroot__ its, itroot__ ite)
Definition: polynomial.hpp:165
data_type coefficient(const index_type &i) const
Definition: polynomial.hpp:77
void adjust_zero(const data_type &small)
Definition: polynomial.hpp:108
Definition: math.hpp:20
polynomial< data_type > * fppp() const
Definition: polynomial.hpp:324
Derived1__::Scalar distance(const Eigen::MatrixBase< Derived1__ > &p1, const Eigen::MatrixBase< Derived2__ > &p2)
Definition: distance.hpp:33
polynomial< data_type > * fp() const
Definition: polynomial.hpp:298
polynomial(const coefficient_type &c)
Definition: polynomial.hpp:43
polynomial(const polynomial< data__ > &p)
Definition: polynomial.hpp:44
polynomial< data_type > * f() const
Definition: polynomial.hpp:290
void set_coefficients(const coefficient_type &ain)
Definition: polynomial.hpp:81
polynomial(const data_type &root1, const data_type &root2, const data_type &root3, const data_type &root4)
Definition: polynomial.hpp:48
polynomial(const data_type &root1, const data_type &root2, const data_type &root3)
Definition: polynomial.hpp:47
void set_roots(const data_type &root1, const data_type &root2, const data_type &root3)
Definition: polynomial.hpp:137
void set_roots(const data_type &root1, const data_type &root2, const data_type &root3, const data_type &root4)
Definition: polynomial.hpp:148
data_type fp(const data_type &t) const
Definition: polynomial.hpp:221
polynomial(const data_type &root1, const data_type &root2)
Definition: polynomial.hpp:46
coefficient_type a
Definition: polynomial.hpp:39
void multiply(const polynomial< data_type > &p1, const polynomial< data_type > &p2)
Definition: polynomial.hpp:383
data_type f(const data_type &t) const
Definition: polynomial.hpp:199
Definition: polynomial.hpp:31
polynomial(const data_type &root)
Definition: polynomial.hpp:45
Eigen::Matrix< data_type, Eigen::Dynamic, 1 > coefficient_type
Definition: polynomial.hpp:35
polynomial< data__ > & operator=(const polynomial< data__ > &p)
Definition: polynomial.hpp:62
void set_roots(const data_type &root)
Definition: polynomial.hpp:118
polynomial< data__ > & operator=(const data_type &d)
Definition: polynomial.hpp:52
index_type degree() const
Definition: polynomial.hpp:72
void subtract(const polynomial< data_type > &p1, const polynomial< data_type > &p2)
Definition: polynomial.hpp:360
void multiply(const data_type &d)
Definition: polynomial.hpp:404
void add(const polynomial< data_type > &p1, const polynomial< data_type > &p2)
Definition: polynomial.hpp:337
polynomial(itroot__ its, itroot__ ite)
Definition: polynomial.hpp:50
void negative()
Definition: polynomial.hpp:481
data_type fpp(const data_type &t) const
Definition: polynomial.hpp:244
index_type size() const
Definition: polynomial.hpp:488
void get_coefficients(coefficient_type &aout) const
Definition: polynomial.hpp:86
polynomial()
Definition: polynomial.hpp:42
void divide(polynomial< data_type > &prem, const polynomial< data_type > &p1, const polynomial< data_type > &p2)
Definition: polynomial.hpp:410
data_type fppp(const data_type &t) const
Definition: polynomial.hpp:267
coefficient_type::Index index_type
Definition: polynomial.hpp:36
void divide(const data_type &d)
Definition: polynomial.hpp:475
bool next_combination(const it__ itb, it__ itk, const it__ ite, comp__ comp)
Definition: combination.hpp:35
void set_roots(const data_type &root1, const data_type &root2)
Definition: polynomial.hpp:127
void compress()
Definition: polynomial.hpp:91
data__ data_type
Definition: polynomial.hpp:34
polynomial< data_type > * fpp() const
Definition: polynomial.hpp:311