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