OpenWalnut  1.3.1
WPolynomialEquationSolvers.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WPOLYNOMIALEQUATIONSOLVERS_H
26 #define WPOLYNOMIALEQUATIONSOLVERS_H
27 
28 #include <complex>
29 #include <iostream>
30 #include <sstream>
31 #include <utility>
32 
33 #include "../exceptions/WEquationHasNoRoots.h"
34 #include "../WLimits.h"
35 
36 
37 /**
38  * Solves a quadratic equation: ax^2 + bx + c = 0 where a,b,c are real coefficient.
39  *
40  * \param a first coefficient
41  * \param b second coefficient
42  * \param c remainder
43  *
44  * \throw WEquationHasNoRoots In case there are no roots for this equation.
45  *
46  * \return Solutions to the equation above as complex numbers. If only one solution exists both
47  * complex numbers are the same!
48  */
49 template< typename T > typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c );
50 
51 
52 
53 template< typename T >
54 inline typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c )
55 {
56  typename std::pair< typename std::complex< T >, typename std::complex< T > > result( std::complex< T >( 0.0, 0.0 ), std::complex< T >( 0.0, 0.0 ) ); // NOLINT line length
57  if( a != 1.0 ) // normalize
58  {
59  if( std::abs( a ) <= wlimits::DBL_EPS ) // a ≃ 0.0
60  {
61  if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0
62  {
63  if( std::abs( c ) > wlimits::DBL_EPS ) // there is no solution
64  {
65  std::stringstream ss;
66  ss << "The equation: " << a << "x^2 + " << b << "x + " << c << " = 0.0 has no solutions!";
67  throw WEquationHasNoRoots( ss.str() );
68  }
69  else // all coefficents are zero so we return 0.0 as a solution
70  {
71  return result;
72  }
73  }
74  else // Note equation degrades to linear equation:
75  {
76  result.first = std::complex< T >( -c / b, 0.0 );
77  result.second = result.first;
78  return result;
79  }
80  }
81  else
82  {
83  b /= a;
84  c /= a;
85  a = 1.0;
86  }
87  }
88 
89  // equation is in normal form
90  double discriminant = b * b - 4.0 * c;
91  if( discriminant < 0.0 )
92  {
93  result.first = std::complex< T >( -b / 2.0, std::sqrt( std::abs( discriminant ) ) / 2.0 );
94  result.second = std::complex< T >( -b / 2.0, -std::sqrt( std::abs( discriminant ) ) / 2.0 );
95  }
96  else if( discriminant > 0.0 )
97  {
98  result.first = std::complex< T >( -b / 2.0 + std::sqrt( discriminant ) / 2.0, 0.0 );
99  result.second = std::complex< T >( -b / 2.0 - std::sqrt( discriminant ) / 2.0 , 0.0 );
100  }
101  else // discriminant ≃ 0.0
102  {
103  // NOTE: if b==-0.0 and discriminant==0.0 then result.first and result.second would be different. To prevent this we check if b ≃ 0.0
104  if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0
105  {
106  result.first = std::complex< T >( 0.0, 0.0 );
107  result.second = result.first;
108  }
109  else
110  {
111  result.first = std::complex< T >( -b / 2, 0.0 );
112  result.second = result.first;
113  }
114  }
115  return result;
116 }
117 
118 #endif // WPOLYNOMIALEQUATIONSOLVERS_H