]>
Commit | Line | Data |
---|---|---|
1e59de90 TL |
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 <algorithm> | |
8 | #include <array> | |
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 | |
16 | // complex roots. Follows Numerical Recipes, Chapter 5, section 6. NB: A better | |
17 | // algorithm apparently exists: Algorithm 954: An Accurate and Efficient Cubic | |
18 | // and Quartic Equation Solver for Physical Applications However, I don't have | |
19 | // access to that paper! | |
20 | template <typename Real> | |
21 | std::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 we | |
76 | // only have one real root if R^2 >= Q^3. But this isn't true; we can | |
77 | // even see this from equation 5.6.18. The condition for having three | |
78 | // real roots is that A = B. It *is* the case that if we're in this | |
79 | // branch, and we have 3 real roots, two are a double root. Take | |
80 | // (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double | |
81 | // 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 the | |
91 | // 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 single | |
102 | // 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 residual | |
123 | // eps*|rp'(r)| (second element) for a root. Recall that for a numerically | |
124 | // computed root r satisfying r = r_0(1+eps) of a function p, |p(r)| <= | |
125 | // eps|rp'(r)|. | |
126 | template <typename Real> | |
127 | std::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 are | |
141 | // 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, A | |
152 | // Graduate Introduction to Numerical Methods, Section 3.2.1. | |
153 | template <typename Real> | |
154 | Real 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::tools | |
176 | #endif |