Finish up cubic roots.

doc/roots/cubic_roots.qbk

@@ -0,0 +1,78 @@
+[/
+Copyright (c) 2021 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:cubic_roots Roots of Cubic Polynomials]
+
+[heading Synopsis]
+
+```
+#include <boost/math/roots/cubic_roots.hpp>
+
+namespace boost::math::tools {
+
+// Solves ax³ + bx² + cx + d = 0.
+std::array<Real,3> cubic_roots(Real a, Real b, Real c, Real d);
+}
+```
+
+[heading Background]
+
+The `cubic_roots` function extracts all real roots of a cubic polynomial ax³ + bx² + cx + d.
+The result is a `std::array<Real, 3>`, which has length three, irrespective of whether there are 3 real roots.
+There is always 1 real root, and hence (barring overflow or other exceptional circumstances), the first element of the
+`std::array` is always populated.
+If there is only one real root of the polynomial, then the second and third elements are set to `nans`.
+The roots are returned in nondecreasing order.
+
+Be careful with double roots.
+First, if you have a real double root, it is numerically indistinguishable from a complex conjugate pair of roots,
+where the complex part is tiny.
+Second, the condition number of rootfinding is infinite at a double root,
+so even changes as subtle as different instruction generation can change the outcome.
+We have some heuristics in place to detect double roots, but these should be regarded with suspicion.
+
+
+[heading Example]
+
+```
+#include <iostream>
+#include <sstream>
+#include <boost/math/tools/cubic_roots.hpp>
+
+using boost::math::tools::cubic_roots;
+
+template<typename Real>
+std::string print_roots(std::array<Real, 3> const & roots) {
+    std::ostringstream out;
+    out << "{" << roots[0] << ", " << roots[1] << ", " << roots[2] << "}";
+    return out.str();
+}
+
+int main() {
+    // Solves x³ - 6x² + 11x - 6 = (x-1)(x-2)(x-3).
+    auto roots = cubic_roots(1.0, -6.0, 11.0, -6.0);
+    std::cout << "The roots of  x³ - 6x² + 11x - 6 are " << print_roots(roots) << ".\n";
+
+    // Double root: YMMV:
+    // (x+1)²(x-2) = x³ - 3x - 2:
+    roots = cubic_roots(1.0, 0.0, -3.0, -2.0);
+    std::cout << "The roots of  x³ - 3x - 2 are " << print_roots(roots) << ".\n";
+}
+```
+
+[heading Performance and Accuracy]
+
+On an Intel laptop chip running at 2700 MHz, we observe 3 roots taking ~90ns to compute.
+If the polynomial only possesses a single real root, it takes ~50ns.
+See `reporting/performance/cubic_roots_performance.cpp`.
+
+The forward error cannot be effectively bounded.
+However, the residuals should be constrained by p(x) <= ε|xp'(x)|.
+
+
+[endsect]
+[/section:cubic_roots]

doc/roots/roots_overview.qbk

@@ -17,6 +17,7 @@ There are several fully-worked __root_finding_examples, including:
 
 [include roots_without_derivatives.qbk]
 [include roots.qbk]
+[include cubic_roots.qbk]
 [include root_finding_examples.qbk]
 [include minima.qbk]
 [include root_comparison.qbk]

include/boost/math/tools/cubic_roots.hpp

@@ -73,45 +73,33 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
         roots[0] = -2*rtQ*ct - p/3;
         roots[1] = -rtQ*(-ct + sqrt(Real(3))*st) - p/3;
         roots[2] = rtQ*(ct + sqrt(Real(3))*st) - p/3;
-        // This formula is not super accurate.
-        // Do a cleanup iteration.
-        for (auto & r : roots) {
-            // Horner's method.
-            // Here I'll take John Gustaffson's opinion that the fma is a *distinct* operation from a*x +b:
-            // Make sure to compile these fmas into a single instruction!
-            Real f = fma(a, r, b);
-            f = fma(f,r,c);
-            f = fma(f,r,d);
-            Real df = fma(3*a, r, 2*b);
-            df = fma(df, r, c);
-            if (df != 0) {
-                // No standard library feature for fused-divide add!
-                r -= f/df;
-            }
+    } else {
+        // In Numerical Recipes, Chapter 5, Section 6, it is claimed that we only have one real root
+        // if R^2 >= Q^3. But this isn't true; we can even see this from equation 5.6.18.
+        // The condition for having three real roots is that A = B.
+        // It *is* the case that if we're in this branch, and we have 3 real roots, two are a double root.
+        // Take (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double root at x = -1,
+        // and it gets sent into this branch.
+        Real arg = R*R - Q*Q*Q;
+        Real A = -detail::sgn(R)*cbrt(abs(R) + sqrt(arg));
+        Real B = 0;
+        if (A != 0) {
+            B = Q/A;
+        }
+        roots[0] = A + B - p/3;
+        // Yes, we're comparing floats for equality:
+        // Any perturbation pushes the roots into the complex plane; out of the bailiwick of this routine.
+        if (A == B || arg == 0) {
+            roots[1] = -A - p/3;
+            roots[2] = -A - p/3;
         }
-        std::sort(roots.begin(), roots.end());
-        return roots;
-    }
-    // In Numerical Recipes, Chapter 5, Section 6, it is claimed that we only have one real root
-    // if R^2 >= Q^3. But this isn't true; we can even see this from equation 5.6.18.
-    // The condition for having three real roots is that A = B.
-    // It *is* the case that if we're in this branch, and we have 3 real roots, two are a double root.
-    // Take (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double root at x = -1,
-    // and it gets sent into this branch.
-    Real arg = R*R - Q*Q*Q;
-    Real A = -detail::sgn(R)*cbrt(abs(R) + sqrt(arg));
-    Real B = 0;
-    if (A != 0) {
-        B = Q/A;
-    }
-    roots[0] = A + B - p/3;
-    // Yes, we're comparing floats for equality:
-    // Any perturbation pushes the roots into the complex plane; out of the bailiwick of this routine.
-    if (A == B || arg == 0) {
-        roots[1] = -A - p/3;
-        roots[2] = -A - p/3;
     }
+    // Root polishing:
     for (auto & r : roots) {
+        // Horner's method.
+        // Here I'll take John Gustaffson's opinion that the fma is a *distinct* operation from a*x +b:
+        // Make sure to compile these fmas into a single instruction and not a function call!
+        // (I'm looking at you Windows.)
         Real f = fma(a, r, b);
         f = fma(f,r,c);
         f = fma(f,r,d);

reporting/performance/cubic_roots_performance.cpp

@@ -16,8 +16,9 @@ template<class Real>
 void CubicRoots(benchmark::State& state)
 {
     std::random_device rd;
-    //auto seed = rd();
-    uint32_t seed = 416683252;
+    auto seed = rd();
+    // This seed generates 3 real roots:
+    //uint32_t seed = 416683252;
     std::mt19937_64 mt(seed);
     std::uniform_real_distribution<Real> unif(-10, 10);
 
@@ -30,12 +31,10 @@ void CubicRoots(benchmark::State& state)
         auto roots = cubic_roots(a,b,c,d);
         benchmark::DoNotOptimize(roots[0]);
     }
-    std::cout << "Just solved " << a << "x^3 + " << b << "x^2 + " << c << "x + " << d << "\n";
-    std::cout << "This was generated by seed " << seed << "\n";
 }
 
-//BENCHMARK_TEMPLATE(CubicRoots, float);
+BENCHMARK_TEMPLATE(CubicRoots, float);
 BENCHMARK_TEMPLATE(CubicRoots, double);
-//BENCHMARK_TEMPLATE(CubicRoots, long double);
+BENCHMARK_TEMPLATE(CubicRoots, long double);
 
 BENCHMARK_MAIN();

test/Jamfile.v2

@@ -990,6 +990,7 @@ test-suite misc :
    [ run signal_statistics_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
    [ run anderson_darling_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
    [ run ljung_box_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
+   [ run cubic_roots_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
    [ run test_t_test.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <define>BOOST_MATH_TEST_FLOAT128 <linkflags>"-Bstatic -lquadmath -Bdynamic" ] [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ]  ]
    [ run test_z_test.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <define>BOOST_MATH_TEST_FLOAT128 <linkflags>"-Bstatic -lquadmath -Bdynamic" ] [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ]  ]
    [ run bivariate_statistics_test.cpp : : : [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ]  ]

test/cubic_roots_test.cpp

@@ -98,13 +98,13 @@ void test_zero_coefficients()
 
         auto roots = cubic_roots(a, b, c, d);
         // I could check the condition number here, but this is fine right?
-        if(!CHECK_ULP_CLOSE(r[0], roots[0], 3)) {
+        if(!CHECK_ULP_CLOSE(r[0], roots[0], 25)) {
             std::cerr << "  Polynomial x^3 + " << b << "x^2 + " << c << "x + " << d << " has roots {";
             std::cerr << r[0] << ", " << r[1] << ", " << r[2] << "}, but the computed roots are {";
             std::cerr << roots[0] << ", " << roots[1] << ", " << roots[2] << "}\n";
         }
-        CHECK_ULP_CLOSE(r[1], roots[1], 3);
-        CHECK_ULP_CLOSE(r[2], roots[2], 3);
+        CHECK_ULP_CLOSE(r[1], roots[1], 25);
+        CHECK_ULP_CLOSE(r[2], roots[2], 25);
     }
 }