polynomial.cpp
Go to the documentation of this file.
1 
12 /*****************************************************************************
13 ** Includes
14 *****************************************************************************/
15 
16 #include <iostream>
17 //#include <ecl/linear_algebra.hpp>
18 #include <ecl/exceptions/standard_exception.hpp>
19 #include <ecl/math/constants.hpp>
20 #include <ecl/math/fuzzy.hpp>
21 #include <ecl/math/simple.hpp> // cube_root
22 #include "../../include/ecl/geometry/polynomial.hpp"
23 
24 /*****************************************************************************
25 ** Namespaces
26 *****************************************************************************/
27 
28 namespace ecl {
29 
30 /*****************************************************************************
31 ** Implementation [Synthetic Division]
32 *****************************************************************************/
33 
34 LinearFunction Division< QuadraticPolynomial >::operator()(const QuadraticPolynomial &p, const double &factor, double &remainder) {
35  double a, b; // linear coefficients
36  a = p.coefficients()[2];
37  b = p.coefficients()[1] + factor*a;
38  LinearFunction f;
39  f.coefficients() << b, a;
40  remainder = p.coefficients()[0]+factor*b;
41  return f;
42 }
43 
44 QuadraticPolynomial Division< CubicPolynomial >::operator()(const CubicPolynomial &p, const double &factor, double &remainder) {
45  double a, b, c; // quadratic coefficients
46  a = p.coefficients()[3];
47  b = p.coefficients()[2] + factor*a;
48  c = p.coefficients()[1] + factor*b;
49  QuadraticPolynomial q;
50  q.coefficients() << c, b, a;
51  remainder = p.coefficients()[0]+factor*c;
52  return q;
53 }
54 
55 /*****************************************************************************
56 ** Implementation [Minimum|Maximum][Polynomial]
57 *****************************************************************************/
58 
59 Array<double> Roots< LinearFunction >::operator()(const LinearFunction& p) {
60  Array<double> intercepts;
61  double a = p.coefficients()[1];
62  double b = p.coefficients()[0];
63  if ( a != 0 ) {
64  intercepts.resize(1);
65  intercepts << -1.0*b/a;
66  }
67  return intercepts;
68 }
69 
70 Array<double> Roots< QuadraticPolynomial >::operator()(const QuadraticPolynomial& p) {
71  Array<double> intercepts;
72  double a = p.coefficients()[2];
73  double b = p.coefficients()[1];
74  double c = p.coefficients()[0];
75  if ( a == 0 ) {
76  LinearFunction f;
77  f.coefficients() << c, b;
78  intercepts = Roots<LinearFunction>()(f);
79  } else {
80  double discriminant = b*b - 4*a*c;
81  if ( discriminant > 0.0 ) {
82  intercepts.resize(2);
83  intercepts << (-b + sqrt(discriminant))/(2*a), (-b - sqrt(discriminant))/(2*a);
84  } else if ( discriminant == 0.0 ) {
85  intercepts.resize(1);
86  intercepts << -b/(2*a);
87  }
88  }
89  return intercepts;
90 }
91 
92  // http://en.wikipedia.org/wiki/Cubic_function#Roots_of_a_cubic_function
93 Array<double> Roots< CubicPolynomial >::operator()(const CubicPolynomial& polynomial) {
94  Array<double> intercepts;
95  double a = polynomial.coefficients()[3];
96  double b = polynomial.coefficients()[2];
97  double c = polynomial.coefficients()[1];
98  double d = polynomial.coefficients()[0];
99 
100 // form the companion matrix (http://stackoverflow.com/questions/2003465/fastest-numerical-solution-of-a-real-cubic-polynomial)
101 // ecl::linear_algebra::Matrix3d A;
102 // A << 0, 0, -d/a, 1, 0, -c/a, 0, 1, -b/a;
103 // ecl::linear_algebra::EigenSolver<ecl::linear_algebra::Matrix3d> eigensolver(A);
104 // if (eigensolver.info() != ecl::linear_algebra::Success) abort();
105 // cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
106 
107  // Monic Trinomial coefficients
108  double p = (3*a*c - b*b)/(3*a*a);
109  double q = (2*b*b*b - 9*a*b*c + 27*a*a*d)/(27*a*a*a);
110  double discriminant = p*p*p/27 + q*q/4;
111 // std::cout << "p: " << p << std::endl;
112 // std::cout << "q: " << q << std::endl;
113  // using transform x = t - b/3a, we can solve t^3+pt+q=0
114  double shift = -b/(3*a);
115  if ( ( p == 0 ) && ( q == 0 ) ) {
116  // single, triple root at t = 0
117  intercepts.resize(1);
118  intercepts << shift;
119 // std::cout << "Single triple root" << std::endl;
120  } else if ( p == 0 ) { // && ( q != 0 )
121  // single root from the cube root function
122 // std::cout << "Single root from cube root function" << std::endl;
123  intercepts.resize(1);
124  intercepts << ecl::cube_root(-q)+ shift;
125  } else if ( q == 0 ) {
126 // std::cout << "Three real roots" << std::endl;
127  // three real roots
128  intercepts.resize(3);
129  intercepts << shift, sqrt(-1*p) + shift, -sqrt(-1*p) + shift;
130  } else if ( discriminant == 0 ) { // && ( p != 0 ) ) {
131 // std::cout << "Double root" << std::endl;
132 // // double root and simple root
133  intercepts.resize(2);
134  intercepts << 3*q/p + shift, (-3*q)/(2*p) + shift;
135  } else if ( discriminant >= 0 ) {
136 // std::cout << "Discriminant: " << discriminant << std::endl;
137  double u = ecl::cube_root(-q/2 + sqrt(discriminant));
138  double v = ecl::cube_root(-q/2 - sqrt(discriminant));
139  intercepts.resize(1);
140  intercepts << u + v + shift;
141  } else { // discriminant < 0 this is cardano's casus irreducibilis and there is three real roots
142  // switch to transcendental solutions (only works for three real roots)
143  // http://en.wikipedia.org/wiki/Cubic_function#Trigonometric_.28and_hyperbolic.29_method
144 // std::cout << "Discriminant: " << discriminant << std::endl;
145  double t_1 = 2.0*sqrt(-p/3.0)*cos((1.0/3.0)*acos( ((3.0*q)/(2.0*p)) * sqrt(-3.0/p)));
146  double t_2 = 2.0*sqrt(-p/3.0)*cos((1.0/3.0)*acos(( (3.0*q)/(2.0*p)) * sqrt(-3.0/p))-(2.0*ecl::pi)/3.0);
147  double t_3 = 2.0*sqrt(-p/3.0)*cos((1.0/3.0)*acos(( (3.0*q)/(2.0*p)) * sqrt(-3.0/p))-(4.0*ecl::pi)/3.0);
148  intercepts.resize(3);
149  intercepts << t_1+shift, t_2+shift, t_3+shift;
150  }
151  return intercepts;
152 }
153 
154 CartesianPoint2d Intersection< LinearFunction >::operator()(const LinearFunction &f, const LinearFunction &g) {
155  CartesianPoint2d point;
156  double a_0 = f.coefficients()[0];
157  double b_0 = f.coefficients()[1];
158  double a_1 = g.coefficients()[0];
159  double b_1 = g.coefficients()[1];
160  if ( isApprox(b_0, b_1) ) { // should use some epsilon distance here.
161  last_operation_failed = true;
162  ecl_throw(StandardException(LOC,OutOfRangeError,"Functions are collinear, no intersection possible."));
163  } else {
164  point.x((a_0 - a_1)/(b_1 - b_0));
165  point.y(f(point.x()));
166  }
167  return point;
168 }
169 
170 double Maximum< LinearFunction >::operator()(
171  const double& x_begin, const double& x_end, const LinearFunction &function) {
172  double max = function(x_begin);
173  double test_max = function(x_end);
174  if ( test_max > max ) {
175  max = test_max;
176  }
177  return max;
178 }
179 
180 double Maximum<CubicPolynomial>::operator()(
181  const double& x_begin, const double& x_end, const CubicPolynomial& cubic) {
182  // 3a_3x^2 + 2a_2x + a_1 = 0
183  double max = cubic(x_begin);
184  double test_max = cubic(x_end);
185  if ( test_max > max ) {
186  max = test_max;
187  }
188  CubicPolynomial::Coefficients coefficients = cubic.coefficients();
189  double a = 3*coefficients[3];
190  double b = 2*coefficients[2];
191  double c = coefficients[1];
192  if ( a == 0 ) {
193  double root = -c/b;
194  if ( ( root > x_begin ) && ( root < x_end ) ) {
195  test_max = cubic(root);
196  if ( test_max > max ) {
197  max = test_max;
198  }
199  }
200  } else {
201  double sqrt_term = b*b-4*a*c;
202  if ( sqrt_term > 0 ) {
203  double root = ( -b + sqrt(b*b-4*a*c))/(2*a);
204  if ( ( root > x_begin ) && ( root < x_end ) ) {
205  test_max = cubic(root);
206  if ( test_max > max ) {
207  max = test_max;
208  }
209  }
210  root = ( -b - sqrt(b*b-4*a*c))/(2*a);
211  if ( ( root > x_begin ) && ( root < x_end ) ) {
212  test_max = cubic(root);
213  if ( test_max > max ) {
214  max = test_max;
215  }
216  }
217  }
218  }
219  return max;
220 }
221 
222 double Minimum< LinearFunction >::operator()(
223  const double& x_begin, const double& x_end, const LinearFunction &function) {
224  double min = function(x_begin);
225  double test_min = function(x_end);
226  if ( test_min < min ) {
227  min = test_min;
228  }
229  return min;
230 }
231 
232 double Minimum<CubicPolynomial>::operator()(
233  const double& x_begin, const double& x_end, const CubicPolynomial& cubic) {
234  // 3a_3x^2 + 2a_2x + a_1 = 0
235  double min = cubic(x_begin);
236  double test_min = cubic(x_end);
237  if ( test_min < min ) {
238  min = test_min;
239  }
240  CubicPolynomial::Coefficients coefficients = cubic.coefficients();
241  double a = 3*coefficients[3];
242  double b = 2*coefficients[2];
243  double c = coefficients[1];
244  if ( a == 0 ) {
245  double root = -c/b;
246  if ( ( root > x_begin ) && ( root < x_end ) ) {
247  test_min = cubic(root);
248  if ( test_min < min ) {
249  min = test_min;
250  }
251  }
252  } else {
253  double sqrt_term = b*b-4*a*c;
254  if ( sqrt_term > 0 ) {
255  double root = ( -b + sqrt(b*b-4*a*c))/(2*a);
256  if ( ( root > x_begin ) && ( root < x_end ) ) {
257  test_min = cubic(root);
258  if ( test_min < min ) {
259  min = test_min;
260  }
261  }
262  root = ( -b - sqrt(b*b-4*a*c))/(2*a);
263  if ( ( root > x_begin ) && ( root < x_end ) ) {
264  test_min = cubic(root);
265  if ( test_min < min ) {
266  min = test_min;
267  }
268  }
269  }
270  }
271  return min;
272 }
273 
274 } // namespace ecl
ecl::Maximum
Primary template functor for the maximum of a continuous function.
Definition: function_math.hpp:60
ecl
Embedded control libraries.
ecl_throw
#define ecl_throw(exception)
ecl::Roots
Primary template functor for the roots of a function (x-axis intercepts).
Definition: function_math.hpp:120
LOC
#define LOC
ecl::pi
double const pi
ecl::Roots< LinearFunction >
X axis intercepts for linear functions.
Definition: polynomial.hpp:1167
ecl::Division
Primary template functor for polynomial division.
Definition: function_math.hpp:107
ecl::OutOfRangeError
OutOfRangeError
ecl::Intersection
Primary template functor for the intersection of like functions.
Definition: function_math.hpp:94
ecl::Polynomial::coefficients
Coefficients & coefficients()
Handle to the coefficient array, use to initialise the polynomial.
Definition: polynomial.hpp:233
f
void f(int i)
ecl::cube_root
Scalar cube_root(const Scalar &x)
ecl::isApprox
bool isApprox(const Scalar &x, const OtherScalar &y, typename numeric_limits< Scalar >::Precision precision=numeric_limits< Scalar >::dummy_precision)
ecl::CartesianPoint
Generic container storing a cartesian point of dimension N.
Definition: cartesian_point.hpp:46
standard_exception.hpp
ecl::Polynomial
Representation of a polynomial function of n-th degree.
Definition: polynomial.hpp:44
ecl::Minimum
Primary template functor for the minimum of a continuous function.
Definition: function_math.hpp:80

ecl_geometry
Author(s): Daniel Stonier
autogenerated on Mon Feb 28 2022 22:18:49