OpenWalnut 1.2.5
|
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 <iostream> 00030 #include <sstream> 00031 #include <utility> 00032 00033 #include "../exceptions/WEquationHasNoRoots.h" 00034 #include "../WLimits.h" 00035 00036 00037 /** 00038 * Solves a quadratic equation: ax^2 + bx + c = 0 where a,b,c are real coefficient. 00039 * 00040 * \param a first coefficient 00041 * \param b second coefficient 00042 * \param c remainder 00043 * 00044 * \throw WEquationHasNoRoots In case there are no roots for this equation. 00045 * 00046 * \return Solutions to the equation above as complex numbers. If only one solution exists both 00047 * complex numbers are the same! 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 ) ); // NOLINT line length 00057 if( a != 1.0 ) // normalize 00058 { 00059 if( std::abs( a ) <= wlimits::DBL_EPS ) // a ≃ 0.0 00060 { 00061 if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0 00062 { 00063 if( std::abs( c ) > wlimits::DBL_EPS ) // there is no solution 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 // all coefficents are zero so we return 0.0 as a solution 00070 { 00071 return result; 00072 } 00073 } 00074 else // Note equation degrades to linear equation: 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 // equation is in normal form 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 // discriminant ≃ 0.0 00102 { 00103 // 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 00104 if( std::abs( b ) <= wlimits::DBL_EPS ) // b ≃ 0.0 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