ceph/src/boost/boost/math/tools/cubic_roots.hpp

   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