OpenWalnut
1.4.0
|
00001 //--------------------------------------------------------------------------- 00002 // 00003 // Project: OpenWalnut ( http://www.openwalnut.org ) 00004 // 00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS 00006 // For more information see http://www.openwalnut.org/copying 00007 // 00008 // This file is part of OpenWalnut. 00009 // 00010 // OpenWalnut is free software: you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as published by 00012 // the Free Software Foundation, either version 3 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // OpenWalnut is distributed in the hope that it will be useful, 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 // GNU Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public License 00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>. 00022 // 00023 //--------------------------------------------------------------------------- 00024 00025 #ifndef WPOLYNOMIALEQUATIONSOLVERS_H 00026 #define WPOLYNOMIALEQUATIONSOLVERS_H 00027 00028 #include <complex> 00029 #include <sstream> 00030 #include <utility> 00031 00032 #include "../exceptions/WEquationHasNoRoots.h" 00033 #include "../WLimits.h" 00034 00035 00036 /** 00037 * Solves a quadratic equation: ax^2 + bx + c = 0 where a,b,c are real coefficient. 00038 * 00039 * \param a first coefficient 00040 * \param b second coefficient 00041 * \param c remainder 00042 * 00043 * \throw WEquationHasNoRoots In case there are no roots for this equation. 00044 * 00045 * \return Solutions to the equation above as complex numbers. If only one solution exists both 00046 * complex numbers are the same! 00047 */ 00048 template< typename T > typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c ); 00049 00050 00051 00052 template< typename T > 00053 inline typename std::pair< typename std::complex< T >, typename std::complex< T > > solveRealQuadraticEquation( T a, T b, T c ) 00054 { 00055 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 00056 if( a != 1.0 ) // normalize 00057 { 00058 if( std::abs( a ) <= wlimits::DBL_EPS ) // a ≃ 0.0 00059 { 00060 if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0 00061 { 00062 if( std::abs( c ) > wlimits::DBL_EPS ) // there is no solution 00063 { 00064 std::stringstream ss; 00065 ss << "The equation: " << a << "x^2 + " << b << "x + " << c << " = 0.0 has no solutions!"; 00066 throw WEquationHasNoRoots( ss.str() ); 00067 } 00068 else // all coefficents are zero so we return 0.0 as a solution 00069 { 00070 return result; 00071 } 00072 } 00073 else // Note equation degrades to linear equation: 00074 { 00075 result.first = std::complex< T >( -c / b, 0.0 ); 00076 result.second = result.first; 00077 return result; 00078 } 00079 } 00080 else 00081 { 00082 b /= a; 00083 c /= a; 00084 a = 1.0; 00085 } 00086 } 00087 00088 // equation is in normal form 00089 double discriminant = b * b - 4.0 * c; 00090 if( discriminant < 0.0 ) 00091 { 00092 result.first = std::complex< T >( -b / 2.0, std::sqrt( std::abs( discriminant ) ) / 2.0 ); 00093 result.second = std::complex< T >( -b / 2.0, -std::sqrt( std::abs( discriminant ) ) / 2.0 ); 00094 } 00095 else if( discriminant > 0.0 ) 00096 { 00097 result.first = std::complex< T >( -b / 2.0 + std::sqrt( discriminant ) / 2.0, 0.0 ); 00098 result.second = std::complex< T >( -b / 2.0 - std::sqrt( discriminant ) / 2.0 , 0.0 ); 00099 } 00100 else // discriminant ≃ 0.0 00101 { 00102 // 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 00103 if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0 00104 { 00105 result.first = std::complex< T >( 0.0, 0.0 ); 00106 result.second = result.first; 00107 } 00108 else 00109 { 00110 result.first = std::complex< T >( -b / 2, 0.0 ); 00111 result.second = result.first; 00112 } 00113 } 00114 return result; 00115 } 00116 00117 #endif // WPOLYNOMIALEQUATIONSOLVERS_H