Code-Eli  0.3.6
sturm_count.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_root_sturm_count_hpp
14 #define eli_mutil_poly_root_sturm_count_hpp
15 
16 #include <cmath>
17 #include <vector>
18 
19 #include "eli/code_eli.hpp"
20 
24 
25 namespace eli
26 {
27  namespace mutil
28  {
29  namespace poly
30  {
31  namespace root
32  {
33  template<typename data__>
34  void sturm_functions(std::vector< polynomial<data__> > &sturm_fun, const polynomial<data__> &f)
35  {
36  size_t i, n;
37 
38  // resize the Sturm function vector
39  n=f.degree()+1;
40  sturm_fun.resize(n);
41 
42  polynomial<data__> *ptemp(f.fp()), pjunk;
43 
44  // initialize the Sturm functions
45  sturm_fun[0]=f;
46  ptemp=f.fp();
47  sturm_fun[1]=*ptemp;
48  delete ptemp;
49 
50  // build the Sturm functions
51  sturm_fun[0].divide(std::abs(sturm_fun[0].coefficient(sturm_fun[0].degree())));
52  sturm_fun[1].divide(std::abs(sturm_fun[1].coefficient(sturm_fun[1].degree())));
53  for (i=2; i<n; ++i)
54  {
55  if (sturm_fun[i-1].degree()>0)
56  {
57  data__ small_no, max_c(0);
59 
60  sturm_fun[i-1].get_coefficients(c);
61  for (typename polynomial<data__>::index_type j=0; j<=sturm_fun[i-1].degree(); ++j)
62  {
63  if (std::abs(c(j))>max_c)
64  max_c=std::abs(c(j));
65  }
66  small_no=max_c*std::sqrt(std::numeric_limits<data__>::epsilon());
67 
68  pjunk.divide(sturm_fun[i], sturm_fun[i-2], sturm_fun[i-1]);
69  sturm_fun[i].negative();
70  if (std::abs(sturm_fun[i].coefficient(sturm_fun[i].degree()))>small_no)
71  sturm_fun[i].divide(std::abs(sturm_fun[i].coefficient(sturm_fun[i].degree())));
72  sturm_fun[i].adjust_zero(small_no);
73  }
74  else
75  {
76  sturm_fun[i]=0;
77  }
78  }
79  }
80 
81  template<typename data__, typename data2__>
82  int sturm_count(const std::vector< polynomial<data__> > &sturm_fun, const data2__ &xmin, const data2__ &xmax)
83  {
84  std::vector<data__> eval(sturm_fun.size());
85  int count_min, count_max;
86 
87  // evaluate functions at end points and count sign changes
88  for (size_t i=0; i<eval.size(); ++i)
89  eval[i]=sturm_fun[i].f(static_cast<data__>(xmin));
90  count_min = eli::mutil::poly::root::sign_changes(eval.begin(), eval.end());
91  for (size_t i=0; i<eval.size(); ++i)
92  eval[i]=sturm_fun[i].f(static_cast<data__>(xmax));
93  count_max = eli::mutil::poly::root::sign_changes(eval.begin(), eval.end());
94 
95  // return the difference between sign changes for min value and sign changes for max value
96  return count_min-count_max;
97  }
98 
99  template<typename data__, typename data2__>
100  int sturm_count(const polynomial<data__> &f, const data2__ &xmin, const data2__ &xmax)
101  {
102  // short circuit degenerate cases
103  if (xmax<=xmin)
104  return 0;
105 
106  // catch another degenerate case
107  if (f.degree()==0)
108  return 0;
109 
110  std::vector< polynomial<data__> > sturm_fun;
111 
112  // calculate the Sturm functions
113  sturm_functions(sturm_fun, f);
114 
115  return sturm_count(sturm_fun, xmin, xmax);
116  }
117 
118  template<typename data__>
119  int sturm_count(const polynomial<data__> &f, bool positive)
120  {
121  typename polynomial<data__>::data_type xmin, xmax;
122 
123  // find appropriate xmin and xmax
124  if (positive)
125  {
126  xmin=static_cast<typename polynomial<data__>::data_type>(0);
128  }
129  else
130  {
132  xmax=static_cast<typename polynomial<data__>::data_type>(0);
133  }
134 
135  return sturm_count(f, xmin, xmax);
136  }
137 
138  template<typename data__>
139  int sturm_count(const std::vector< polynomial<data__> > &sturm_fun, bool positive)
140  {
141  typename polynomial<data__>::data_type xmin, xmax;
142 
143  // find appropriate xmin and xmax
144  if (positive)
145  {
146  xmin=static_cast<typename polynomial<data__>::data_type>(0);
147  xmax=eli::mutil::poly::root::max_radius(sturm_fun[0]);
148  }
149  else
150  {
151  xmin=-eli::mutil::poly::root::max_radius(sturm_fun[0]);
152  xmax=static_cast<typename polynomial<data__>::data_type>(0);
153  }
154 
155  return sturm_count(sturm_fun, xmin, xmax);
156  }
157 
158  template<typename data__>
160  {
161  typename polynomial<data__>::data_type xmin, xmax;
162 
163  // find appropriate xmin and xmax
165  xmin=-xmax;
166 
167  return sturm_count(f, xmin, xmax);
168  }
169 
170  template<typename data__>
171  int sturm_count(const std::vector< polynomial<data__> > &sturm_fun)
172  {
173  typename polynomial<data__>::data_type xmin, xmax;
174 
175  // find appropriate xmin and xmax
176  xmax=eli::mutil::poly::root::max_radius(sturm_fun[0]);
177  xmin=-xmax;
178 
179  return sturm_count(sturm_fun, xmin, xmax);
180  }
181  }
182  }
183  }
184 }
185 
186 #endif
data__ max_radius(const polynomial< data__ > &f)
Definition: radius.hpp:31
void adjust_zero(const data_type &small)
Definition: polynomial.hpp:108
Definition: math.hpp:20
data_type fp(const data_type &t) const
Definition: polynomial.hpp:221
Definition: polynomial.hpp:31
Eigen::Matrix< data_type, Eigen::Dynamic, 1 > coefficient_type
Definition: polynomial.hpp:35
int sturm_count(const std::vector< polynomial< data__ > > &sturm_fun, const data2__ &xmin, const data2__ &xmax)
Definition: sturm_count.hpp:82
index_type degree() const
Definition: polynomial.hpp:72
void negative()
Definition: polynomial.hpp:481
void sturm_functions(std::vector< polynomial< data__ > > &sturm_fun, const polynomial< data__ > &f)
Definition: sturm_count.hpp:34
int sign_changes(const it__ itb, const it__ ite)
Definition: sign_changes.hpp:28
void get_coefficients(coefficient_type &aout) const
Definition: polynomial.hpp:86
void divide(polynomial< data_type > &prem, const polynomial< data_type > &p1, const polynomial< data_type > &p2)
Definition: polynomial.hpp:410
coefficient_type::Index index_type
Definition: polynomial.hpp:36
data__ data_type
Definition: polynomial.hpp:34