OpenWalnut 1.2.5
WPolynomialEquationSolvers.h
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
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends