roots.hpp

Go to the documentation of this file.

/** \file roots.hpp
 * \author Jared R. Males (jaredmales@gmail.com)
 * \brief Declares and defines functions for finding roots
 * \ingroup gen_math_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// mxlib is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef math_roots_hpp
#define math_roots_hpp
 
#include <vector>
#include <complex>
#include <cmath>
 
namespace mx
{
namespace math
{
 
/// Calculate the descriminant for the general cubic
/**
 * \returns the descriminant
 */
template <typename realT>
realT cubicDescriminant( const realT &a, ///< [in] the 3rd order coefficient of the general cubic
                         const realT &b, ///< [in] the 2nd order coefficient of the general cubic
                         const realT &c, ///< [in] the 1st order coefficient of the general cubic
                         const realT &d  ///< [in] the 0th order coefficient of the general cubic
)
{
    return 18 * a * b * c * d - 4 * b * b * b * d + b * b * c * c - 4 * a * c * c * c - 27 * a * a * d * d;
}
 
/// Calculate the descriminant for the depressed cubic
/**
 * \returns the descriminant
 */
template <typename realT>
realT cubicDescriminant( const realT &p, ///< [in] the 0th order coefficient of the depressed cubic
                         const realT &q  ///< [in] the 0th order coefficient of the depressed cubic
)
{
    return -1 * ( 4 * p * p * p + 27 * q * q );
}
 
/// Convert a general cubic equation to depressed form
/** The general cubic is
 * \f[
 * ax^3 + bx^2 + cx + d = 0
 * \f]
 * which can be converted to compressed form
 * \f[
 * t^3 + pt + q = 0
 * \f]
 *
 */
template <typename realT>
void cubicDepressed( realT &p,       ///< [out] the 1st order coefficient of the depressed cubic
                     realT &q,       ///< [out] the 0th order coefficient of the depressed cubic
                     const realT &a, ///< [in] the 3rd order coefficient of the general cubic
                     const realT &b, ///< [in] the 2nd order coefficient of the general cubic
                     const realT &c, ///< [in] the 1st order coefficient of the general cubic
                     const realT &d  ///< [in] the 0th order coefficient of the general cubic
)
{
    p = ( 3 * a * c - b * b ) / ( 3 * a * a );
    q = ( 2 * b * b * b - 9 * a * b * c + 27 * a * a * d ) / ( 27 * a * a * a );
}
 
/// Calculate the real root for a depressed cubic with negative descriminant
/**
 * \returns the real root for a depressed cubic defined by \p p and \p q
 *
 * \throws std::runtime_error if there are 3 real roots
 *
 */
template <typename realT>
realT cubicRealRoot( const realT &p, ///< [in] the 1st order coefficient of the depressed cubic
                     const realT &q  ///< [in] the 0th order coefficient of the depressed cubic
)
{
    realT D = q * q / 4 + p * p * p / 27;
 
    if( D < 0 )
        throw std::runtime_error( "cannot apply Cardano's formula, find roots using...." );
 
    D = sqrt( D );
 
    realT u1 = -q / 2 + D;
    realT u2 = -q / 2 - D;
 
    return std::cbrt( u1 ) + std::cbrt( u2 );
}
 
/// Find the roots of the general quartic equation
/** Finds the roots of
  * \f[
    f(x) = a x^4 + b x^3 + c x^2 + d x + e
    \f]
  * using the general formula for quartic roots. See https://en.wikipedia.org/wiki/Quartic_function.
  *
  * \tparam realT is the floating point type used for calculations.
  *
  * \ingroup gen_math
  */
template <typename realT>
void quarticRoots( std::vector<std::complex<realT>> &x, ///< [out] On exit contains the 4 roots, is resized to length 4.
                   realT a,                             ///< [in] the coefficient of the \f$x^4\f$ term.
                   realT b,                             ///< [in] the coefficient of the \f$x^3\f$ term.
                   realT c,                             ///< [in] the coefficient of the \f$x^2\f$ term.
                   realT d,                             ///< [in] the coefficient of the \f$x^1\f$ term.
                   realT e                              ///< [in] the coefficient of the \f$x^0\f$ term.
)
{
    std::complex<realT> p, q, S, Q, Delta0, Delta1;
 
    p = ( 8.0 * a * c - 3.0 * b * b ) / ( 8.0 * a * a );
 
    q = ( b * b * b - 4.0 * a * b * c + 8.0 * a * a * d ) / ( 8.0 * a * a * a );
 
    Delta0 = c * c - 3.0 * b * d + 12.0 * a * e;
    Delta1 = 2.0 * c * c * c - 9.0 * b * c * d + 27.0 * b * b * e + 27.0 * a * d * d - 72.0 * a * c * e;
 
    Q = pow( static_cast<realT>( 0.5 ) *
                 ( Delta1 + sqrt( Delta1 * Delta1 - static_cast<realT>( 4 ) * Delta0 * Delta0 * Delta0 ) ),
             1. / 3. );
 
    S = static_cast<realT>( 0.5 ) *
        sqrt( -static_cast<realT>( 2.0 / 3.0 ) * p + static_cast<realT>( 1.0 / ( 3.0 * a ) ) * ( Q + Delta0 / Q ) );
 
    x.resize( 4 );
 
    x[0] = -b / ( 4 * a ) - S +
           static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p + q / S );
    x[1] = -b / ( 4 * a ) - S -
           static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p + q / S );
    x[2] = -b / ( 4 * a ) + S +
           static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p - q / S );
    x[3] = -b / ( 4 * a ) + S -
           static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p - q / S );
 
} // quarticRoots
 
} // namespace math
} // namespace mx
 
#endif // math_roots_hpp