177 lines · plain
1// (C) Copyright Nick Thompson 2021.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5#ifndef BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP6#define BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP7#include <algorithm>8#include <array>9#include <boost/math/special_functions/sign.hpp>10#include <boost/math/tools/roots.hpp>11 12namespace boost::math::tools {13 14// Solves ax^3 + bx^2 + cx + d = 0.15// Only returns the real roots, as types get weird for real coefficients and16// complex roots. Follows Numerical Recipes, Chapter 5, section 6. NB: A better17// algorithm apparently exists: Algorithm 954: An Accurate and Efficient Cubic18// and Quartic Equation Solver for Physical Applications However, I don't have19// access to that paper!20template <typename Real>21std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {22 using std::abs;23 using std::acos;24 using std::cbrt;25 using std::cos;26 using std::fma;27 using std::sqrt;28 std::array<Real, 3> roots = {std::numeric_limits<Real>::quiet_NaN(),29 std::numeric_limits<Real>::quiet_NaN(),30 std::numeric_limits<Real>::quiet_NaN()};31 if (a == 0) {32 // bx^2 + cx + d = 0:33 if (b == 0) {34 // cx + d = 0:35 if (c == 0) {36 if (d != 0) {37 // No solutions:38 return roots;39 }40 roots[0] = 0;41 roots[1] = 0;42 roots[2] = 0;43 return roots;44 }45 roots[0] = -d / c;46 return roots;47 }48 auto [x0, x1] = quadratic_roots(b, c, d);49 roots[0] = x0;50 roots[1] = x1;51 return roots;52 }53 if (d == 0) {54 auto [x0, x1] = quadratic_roots(a, b, c);55 roots[0] = x0;56 roots[1] = x1;57 roots[2] = 0;58 std::sort(roots.begin(), roots.end());59 return roots;60 }61 Real p = b / a;62 Real q = c / a;63 Real r = d / a;64 Real Q = (p * p - 3 * q) / 9;65 Real R = (2 * p * p * p - 9 * p * q + 27 * r) / 54;66 if (R * R < Q * Q * Q) {67 Real rtQ = sqrt(Q);68 Real theta = acos(R / (Q * rtQ)) / 3;69 Real st = sin(theta);70 Real ct = cos(theta);71 roots[0] = -2 * rtQ * ct - p / 3;72 roots[1] = -rtQ * (-ct + sqrt(Real(3)) * st) - p / 3;73 roots[2] = rtQ * (ct + sqrt(Real(3)) * st) - p / 3;74 } else {75 // In Numerical Recipes, Chapter 5, Section 6, it is claimed that we76 // only have one real root if R^2 >= Q^3. But this isn't true; we can77 // even see this from equation 5.6.18. The condition for having three78 // real roots is that A = B. It *is* the case that if we're in this79 // branch, and we have 3 real roots, two are a double root. Take80 // (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double81 // root at x = -1, and it gets sent into this branch.82 Real arg = R * R - Q * Q * Q;83 Real A = (R >= 0 ? -1 : 1) * cbrt(abs(R) + sqrt(arg));84 Real B = 0;85 if (A != 0) {86 B = Q / A;87 }88 roots[0] = A + B - p / 3;89 // Yes, we're comparing floats for equality:90 // Any perturbation pushes the roots into the complex plane; out of the91 // bailiwick of this routine.92 if (A == B || arg == 0) {93 roots[1] = -A - p / 3;94 roots[2] = -A - p / 3;95 }96 }97 // Root polishing:98 for (auto &r : roots) {99 // Horner's method.100 // Here I'll take John Gustaffson's opinion that the fma is a *distinct*101 // operation from a*x +b: Make sure to compile these fmas into a single102 // instruction and not a function call! (I'm looking at you Windows.)103 Real f = fma(a, r, b);104 f = fma(f, r, c);105 f = fma(f, r, d);106 Real df = fma(3 * a, r, 2 * b);107 df = fma(df, r, c);108 if (df != 0) {109 Real d2f = fma(6 * a, r, 2 * b);110 Real denom = 2 * df * df - f * d2f;111 if (denom != 0) {112 r -= 2 * f * df / denom;113 } else {114 r -= f / df;115 }116 }117 }118 std::sort(roots.begin(), roots.end());119 return roots;120}121 122// Computes the empirical residual p(r) (first element) and expected residual123// eps*|rp'(r)| (second element) for a root. Recall that for a numerically124// computed root r satisfying r = r_0(1+eps) of a function p, |p(r)| <=125// eps|rp'(r)|.126template <typename Real>127std::array<Real, 2> cubic_root_residual(Real a, Real b, Real c, Real d,128 Real root) {129 using std::abs;130 using std::fma;131 std::array<Real, 2> out;132 Real residual = fma(a, root, b);133 residual = fma(residual, root, c);134 residual = fma(residual, root, d);135 136 out[0] = residual;137 138 // The expected residual is:139 // eps*[4|ar^3| + 3|br^2| + 2|cr| + |d|]140 // This can be demonstrated by assuming the coefficients and the root are141 // perturbed according to the rounding model of floating point arithmetic,142 // and then working through the inequalities.143 root = abs(root);144 Real expected_residual = fma(4 * abs(a), root, 3 * abs(b));145 expected_residual = fma(expected_residual, root, 2 * abs(c));146 expected_residual = fma(expected_residual, root, abs(d));147 out[1] = expected_residual * std::numeric_limits<Real>::epsilon();148 return out;149}150 151// Computes the condition number of rootfinding. This is defined in Corless, A152// Graduate Introduction to Numerical Methods, Section 3.2.1.153template <typename Real>154Real cubic_root_condition_number(Real a, Real b, Real c, Real d, Real root) {155 using std::abs;156 using std::fma;157 // There are *absolute* condition numbers that can be defined when r = 0;158 // but they basically reduce to the residual computed above.159 if (root == static_cast<Real>(0)) {160 return std::numeric_limits<Real>::infinity();161 }162 163 Real numerator = fma(abs(a), abs(root), abs(b));164 numerator = fma(numerator, abs(root), abs(c));165 numerator = fma(numerator, abs(root), abs(d));166 Real denominator = fma(3 * a, root, 2 * b);167 denominator = fma(denominator, root, c);168 if (denominator == static_cast<Real>(0)) {169 return std::numeric_limits<Real>::infinity();170 }171 denominator *= root;172 return numerator / abs(denominator);173}174 175} // namespace boost::math::tools176#endif177