solveCubic.h

Go to the documentation of this file.
00001 
00016 #ifndef DLR_NUMERIC_SOLVECUBIC_H
00017 #define DLR_NUMERIC_SOLVECUBIC_H
00018 
00019 #include <complex>
00020 
00021 namespace dlr {
00022 
00023   namespace numeric {
00024     
00052     template <class Type>
00053     bool
00054     solveCubic(Type c0, Type c1, Type c2,
00055                Type& root0, Type& root1, Type& root2);
00056 
00057   
00080     template <class Type>
00081     void
00082     solveCubic(Type c0, Type c1, Type c2,
00083                std::complex<Type>& root0, std::complex<Type>& root1,
00084                std::complex<Type>& root2);
00085 
00086   
00087   } // namespace numeric
00088 
00089 } // namespace dlr
00090 
00091 
00092 /* ======== Inline and template definitions below. ======== */
00093 
00094 #include <cmath>
00095 #include <dlrCommon/constants.h>
00096 
00097 namespace dlr {
00098 
00099   namespace numeric {
00100 
00101     // This function computes the real roots of the cubic polynomial
00102     // x^3 + c0*x^2 + c1*x + c2 = 0.
00103     template <class Type>
00104     bool
00105     solveCubic(Type c0, Type c1, Type c2,
00106                Type& root0, Type& root1, Type& root2)
00107     {
00108       // We follow the formulation in Press et al, "Numerical Recipes,
00109       // The Art of Scientific Computing," third edition, Cambridge
00110       // University Press, 2007.
00111 
00112       bool returnValue = true;
00113       
00114       Type c0Squared = c0 * c0;
00115       Type qq = ((c0Squared - (Type(3.0) * c1)) / Type(9.0));
00116       Type rr = ((Type(2.0) * c0Squared * c0
00117                   - Type(9.0) * c0 * c1
00118                   + Type(27.0) * c2)
00119                  / Type(54.0));
00120 
00121       Type rrSquared = rr * rr;
00122       Type qqCubed = qq * qq * qq;
00123       Type c0OverThree = c0 / Type(3.0);
00124       if(rrSquared < qqCubed) {
00125         // Looks like we have three real roots.
00126         Type theta = std::acos(rr / std::sqrt(qqCubed));
00127         Type minusTwoRootQq = Type(-2.0) * Type(std::sqrt(qq));
00128         Type twoPi = Type(2.0 * dlr::common::constants::pi);
00129 
00130         root0 = (minusTwoRootQq * std::cos(theta / Type(3.0))
00131                  - c0OverThree);
00132         root1 = (minusTwoRootQq * std::cos((theta + twoPi) / Type(3.0))
00133                  - c0OverThree);
00134         root2 = (minusTwoRootQq * std::cos((theta - twoPi) / Type(3.0))
00135                  - c0OverThree);
00136       } else {
00137         // Looks like we have some complex roots.
00138         bool signRr = rr > Type(0.0);
00139         Type absRr = signRr ? rr : -rr;
00140         Type aa = std::pow(absRr + std::sqrt(rrSquared - qqCubed), 1.0 / 3.0);
00141         if(signRr) {
00142           aa = -aa;
00143         }
00144 
00145         Type bb = (aa == Type(0.0)) ? Type(0.0) : (qq / aa);
00146 
00147         root0 = (aa + bb) - c0OverThree;
00148         returnValue = false;
00149       }
00150       return returnValue;
00151     }
00152 
00153 
00154 
00155     // This function computes the (possibly complex) roots of the
00156     // cubic polynomial x^3 + c0*x^2 + c1*x + c2 = 0.
00157     template <class Type>
00158     void
00159     solveCubic(Type c0, Type c1, Type c2,
00160                std::complex<Type>& root0, std::complex<Type>& root1,
00161                std::complex<Type>& root2)
00162     {
00163       // We follow the formulation in Press et al, "Numerical Recipes,
00164       // The Art of Scientific Computing," third edition, Cambridge
00165       // University Press, 2007.
00166       Type c0Squared = c0 * c0;
00167       Type qq = ((c0Squared - (Type(3.0) * c1)) / Type(9.0));
00168       Type rr = ((Type(2.0) * c0Squared * c0
00169                   - Type(9.0) * c0 * c1
00170                   + Type(27.0) * c2)
00171                  / Type(54.0));
00172 
00173       Type rrSquared = rr * rr;
00174       Type qqCubed = qq * qq * qq;
00175       Type c0OverThree = c0 / Type(3.0);
00176       if(rrSquared < qqCubed) {
00177         // Looks like we have three real roots.
00178         Type theta = std::acos(rr / std::sqrt(qqCubed));
00179         Type minusTwoRootQq = Type(-2.0) * Type(std::sqrt(qq));
00180         Type twoPi = Type(2.0 * dlr::common::constants::pi);
00181 
00182         root0.real() = (minusTwoRootQq * std::cos(theta / Type(3.0))
00183                         - c0OverThree);
00184         root0.imag() = 0.0;
00185         root1.real() = (minusTwoRootQq * std::cos((theta + twoPi) / Type(3.0))
00186                         - c0OverThree);
00187         root1.imag() = 0.0;
00188         root2.real() = (minusTwoRootQq * std::cos((theta - twoPi) / Type(3.0))
00189                         - c0OverThree);
00190         root2.imag() = 0.0;
00191       } else {
00192         // Looks like we have some complex roots.
00193         bool signRr = rr > Type(0.0);
00194         Type absRr = signRr ? rr : -rr;
00195         Type aa = std::pow(absRr + std::sqrt(rrSquared - qqCubed), 1.0 / 3.0);
00196         if(signRr) {
00197           aa = -aa;
00198         }
00199 
00200         Type bb = (aa == Type(0.0)) ? Type(0.0) : (qq / aa);
00201 
00202         root0.real() = (aa + bb) - c0OverThree;
00203         root0.imag() = 0.0;
00204         root1.real() = Type(-0.5) * (aa + bb) - c0OverThree;
00205         root1.imag() = (Type(std::sqrt(3.0)) / Type(2.0)) * (aa - bb);
00206         root2.real() = root1.real();
00207         root2.imag() = -root1.imag();
00208       }
00209     }
00210 
00211   } // namespace numeric
00212 
00213 } // namespace dlr
00214 
00215 #endif /* #ifndef DLR_NUMERIC_UTILITIES_H */