1 // (C) Copyright Nick Thompson 2021.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 #ifndef BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP
6 #define BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP
7 #include <array>
8 #include <algorithm>
9 #include <boost/math/special_functions/sign.hpp>
10 #include <boost/math/tools/roots.hpp>
11
12 namespace 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 and complex roots.
16 // Follows Numerical Recipes, Chapter 5, section 6.
17 // NB: A better algorithm apparently exists:
18 // Algorithm 954: An Accurate and Efficient Cubic and Quartic Equation Solver for Physical Applications
19 // However, I don't have access to that paper!
20 template<typename Real>
cubic_roots(Real a,Real b,Real c,Real d)21 std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
22 using std::sqrt;
23 using std::acos;
24 using std::cos;
25 using std::cbrt;
26 using std::abs;
27 using std::fma;
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 we only have one real root
76 // if R^2 >= Q^3. But this isn't true; we can even see this from equation 5.6.18.
77 // The condition for having three real roots is that A = B.
78 // It *is* the case that if we're in this branch, and we have 3 real roots, two are a double root.
79 // Take (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double root at x = -1,
80 // and it gets sent into this branch.
81 Real arg = R*R - Q*Q*Q;
82 Real A = -boost::math::sign(R)*cbrt(abs(R) + sqrt(arg));
83 Real B = 0;
84 if (A != 0) {
85 B = Q/A;
86 }
87 roots[0] = A + B - p/3;
88 // Yes, we're comparing floats for equality:
89 // Any perturbation pushes the roots into the complex plane; out of the bailiwick of this routine.
90 if (A == B || arg == 0) {
91 roots[1] = -A - p/3;
92 roots[2] = -A - p/3;
93 }
94 }
95 // Root polishing:
96 for (auto & r : roots) {
97 // Horner's method.
98 // Here I'll take John Gustaffson's opinion that the fma is a *distinct* operation from a*x +b:
99 // Make sure to compile these fmas into a single instruction and not a function call!
100 // (I'm looking at you Windows.)
101 Real f = fma(a, r, b);
102 f = fma(f,r,c);
103 f = fma(f,r,d);
104 Real df = fma(3*a, r, 2*b);
105 df = fma(df, r, c);
106 if (df != 0) {
107 // No standard library feature for fused-divide add!
108 r -= f/df;
109 }
110 }
111 std::sort(roots.begin(), roots.end());
112 return roots;
113 }
114
115 // Computes the empirical residual p(r) (first element) and expected residual eps*|rp'(r)| (second element) for a root.
116 // Recall that for a numerically computed root r satisfying r = r_0(1+eps) of a function p, |p(r)| <= eps|rp'(r)|.
117 template<typename Real>
cubic_root_residual(Real a,Real b,Real c,Real d,Real root)118 std::array<Real, 2> cubic_root_residual(Real a, Real b, Real c, Real d, Real root) {
119 using std::fma;
120 using std::abs;
121 std::array<Real, 2> out;
122 Real residual = fma(a, root, b);
123 residual = fma(residual,root,c);
124 residual = fma(residual,root,d);
125
126 out[0] = residual;
127
128 Real expected_residual = fma(3*a, root, 2*b);
129 expected_residual = fma(expected_residual, root, c);
130 expected_residual = abs(root*expected_residual)*std::numeric_limits<Real>::epsilon();
131 out[1] = expected_residual;
132 return out;
133 }
134
135 }
136 #endif
137