OpenWalnut  1.4.0
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 <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