00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef WPOLYNOMIALEQUATIONSOLVERS_H
00026 #define WPOLYNOMIALEQUATIONSOLVERS_H
00027
00028 #include <complex>
00029 #include <iostream>
00030 #include <sstream>
00031 #include <utility>
00032
00033 #include "../exceptions/WEquationHasNoRoots.h"
00034 #include "../WLimits.h"
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 template< typename T > typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c );
00050
00051
00052
00053 template< typename T >
00054 inline typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c )
00055 {
00056 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 ) );
00057 if( a != 1.0 )
00058 {
00059 if( std::abs( a ) <= wlimits::DBL_EPS )
00060 {
00061 if( std::abs( b ) <= wlimits::DBL_EPS )
00062 {
00063 if( std::abs( c ) > wlimits::DBL_EPS )
00064 {
00065 std::stringstream ss;
00066 ss << "The equation: " << a << "x^2 + " << b << "x + " << c << " = 0.0 has no solutions!";
00067 throw WEquationHasNoRoots( ss.str() );
00068 }
00069 else
00070 {
00071 return result;
00072 }
00073 }
00074 else
00075 {
00076 result.first = std::complex< T >( -c / b, 0.0 );
00077 result.second = result.first;
00078 return result;
00079 }
00080 }
00081 else
00082 {
00083 b /= a;
00084 c /= a;
00085 a = 1.0;
00086 }
00087 }
00088
00089
00090 double discriminant = b * b - 4.0 * c;
00091 if( discriminant < 0.0 )
00092 {
00093 result.first = std::complex< T >( -b / 2.0, std::sqrt( std::abs( discriminant ) ) / 2.0 );
00094 result.second = std::complex< T >( -b / 2.0, -std::sqrt( std::abs( discriminant ) ) / 2.0 );
00095 }
00096 else if( discriminant > 0.0 )
00097 {
00098 result.first = std::complex< T >( -b / 2.0 + std::sqrt( discriminant ) / 2.0, 0.0 );
00099 result.second = std::complex< T >( -b / 2.0 - std::sqrt( discriminant ) / 2.0 , 0.0 );
00100 }
00101 else
00102 {
00103
00104 if( std::abs( b ) <= wlimits::DBL_EPS )
00105 {
00106 result.first = std::complex< T >( 0.0, 0.0 );
00107 result.second = result.first;
00108 }
00109 else
00110 {
00111 result.first = std::complex< T >( -b / 2, 0.0 );
00112 result.second = result.first;
00113 }
00114 }
00115 return result;
00116 }
00117
00118 #endif // WPOLYNOMIALEQUATIONSOLVERS_H