Line data Source code
1 : /** \file roots.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Declares and defines functions for finding roots
4 : * \ingroup gen_math_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef math_roots_hpp
28 : #define math_roots_hpp
29 :
30 : #include <vector>
31 : #include <complex>
32 : #include <cmath>
33 :
34 : namespace mx
35 : {
36 : namespace math
37 : {
38 :
39 : /// Calculate the descriminant for the general cubic
40 : /**
41 : * \returns the descriminant
42 : */
43 : template <typename realT>
44 : realT cubicDescriminant( const realT &a, ///< [in] the 3rd order coefficient of the general cubic
45 : const realT &b, ///< [in] the 2nd order coefficient of the general cubic
46 : const realT &c, ///< [in] the 1st order coefficient of the general cubic
47 : const realT &d ///< [in] the 0th order coefficient of the general cubic
48 : )
49 : {
50 : 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;
51 : }
52 :
53 : /// Calculate the descriminant for the depressed cubic
54 : /**
55 : * \returns the descriminant
56 : */
57 : template <typename realT>
58 : realT cubicDescriminant( const realT &p, ///< [in] the 0th order coefficient of the depressed cubic
59 : const realT &q ///< [in] the 0th order coefficient of the depressed cubic
60 : )
61 : {
62 : return -1 * ( 4 * p * p * p + 27 * q * q );
63 : }
64 :
65 : /// Convert a general cubic equation to depressed form
66 : /** The general cubic is
67 : * \f[
68 : * ax^3 + bx^2 + cx + d = 0
69 : * \f]
70 : * which can be converted to compressed form
71 : * \f[
72 : * t^3 + pt + q = 0
73 : * \f]
74 : *
75 : */
76 : template <typename realT>
77 0 : void cubicDepressed( realT &p, ///< [out] the 1st order coefficient of the depressed cubic
78 : realT &q, ///< [out] the 0th order coefficient of the depressed cubic
79 : const realT &a, ///< [in] the 3rd order coefficient of the general cubic
80 : const realT &b, ///< [in] the 2nd order coefficient of the general cubic
81 : const realT &c, ///< [in] the 1st order coefficient of the general cubic
82 : const realT &d ///< [in] the 0th order coefficient of the general cubic
83 : )
84 : {
85 0 : p = ( 3 * a * c - b * b ) / ( 3 * a * a );
86 0 : q = ( 2 * b * b * b - 9 * a * b * c + 27 * a * a * d ) / ( 27 * a * a * a );
87 0 : }
88 :
89 : /// Calculate the real root for a depressed cubic with negative descriminant
90 : /**
91 : * \returns the real root for a depressed cubic defined by \p p and \p q
92 : *
93 : * \throws std::runtime_error if there are 3 real roots
94 : *
95 : */
96 : template <typename realT>
97 0 : realT cubicRealRoot( const realT &p, ///< [in] the 1st order coefficient of the depressed cubic
98 : const realT &q ///< [in] the 0th order coefficient of the depressed cubic
99 : )
100 : {
101 0 : realT D = q * q / 4 + p * p * p / 27;
102 :
103 0 : if( D < 0 )
104 0 : throw std::runtime_error( "cannot apply Cardano's formula, find roots using...." );
105 :
106 0 : D = sqrt( D );
107 :
108 0 : realT u1 = -q / 2 + D;
109 0 : realT u2 = -q / 2 - D;
110 :
111 0 : return std::cbrt( u1 ) + std::cbrt( u2 );
112 : }
113 :
114 : /// Find the roots of the general quartic equation
115 : /** Finds the roots of
116 : * \f[
117 : f(x) = a x^4 + b x^3 + c x^2 + d x + e
118 : \f]
119 : * using the general formula for quartic roots. See https://en.wikipedia.org/wiki/Quartic_function.
120 : *
121 : * \tparam realT is the floating point type used for calculations.
122 : *
123 : * \ingroup gen_math
124 : */
125 : template <typename realT>
126 0 : void quarticRoots( std::vector<std::complex<realT>> &x, ///< [out] On exit contains the 4 roots, is resized to length 4.
127 : realT a, ///< [in] the coefficient of the \f$x^4\f$ term.
128 : realT b, ///< [in] the coefficient of the \f$x^3\f$ term.
129 : realT c, ///< [in] the coefficient of the \f$x^2\f$ term.
130 : realT d, ///< [in] the coefficient of the \f$x^1\f$ term.
131 : realT e ///< [in] the coefficient of the \f$x^0\f$ term.
132 : )
133 : {
134 0 : std::complex<realT> p, q, S, Q, Delta0, Delta1;
135 :
136 0 : p = ( 8.0 * a * c - 3.0 * b * b ) / ( 8.0 * a * a );
137 :
138 0 : q = ( b * b * b - 4.0 * a * b * c + 8.0 * a * a * d ) / ( 8.0 * a * a * a );
139 :
140 0 : Delta0 = c * c - 3.0 * b * d + 12.0 * a * e;
141 0 : 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;
142 :
143 0 : Q = pow( static_cast<realT>( 0.5 ) *
144 0 : ( Delta1 + sqrt( Delta1 * Delta1 - static_cast<realT>( 4 ) * Delta0 * Delta0 * Delta0 ) ),
145 0 : 1. / 3. );
146 :
147 0 : S = static_cast<realT>( 0.5 ) *
148 0 : sqrt( -static_cast<realT>( 2.0 / 3.0 ) * p + static_cast<realT>( 1.0 / ( 3.0 * a ) ) * ( Q + Delta0 / Q ) );
149 :
150 0 : x.resize( 4 );
151 :
152 0 : x[0] = -b / ( 4 * a ) - S +
153 0 : static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p + q / S );
154 0 : x[1] = -b / ( 4 * a ) - S -
155 0 : static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p + q / S );
156 0 : x[2] = -b / ( 4 * a ) + S +
157 0 : static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p - q / S );
158 0 : x[3] = -b / ( 4 * a ) + S -
159 0 : static_cast<realT>( 0.5 ) * sqrt( static_cast<realT>( -4 ) * S * S - static_cast<realT>( 2 ) * p - q / S );
160 :
161 0 : } // quarticRoots
162 :
163 : } // namespace math
164 : } // namespace mx
165 :
166 : #endif // math_roots_hpp
|
