math/tools/quartic_roots.hpp

0001 //  (C) Copyright Nick Thompson 2021.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 #ifndef BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP
0006 #define BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP
0007 #include <array>
0008 #include <cmath>
0009 #include <boost/math/tools/cubic_roots.hpp>
0010
0011 namespace boost::math::tools {
0012
0013 namespace detail {
0014
0015 // Make sure the nans are always at the back of the array:
0016 template<typename Real>
0017 bool comparator(Real r1, Real r2) {
0018    using std::isnan;
0019    if (isnan(r1)) { return false; }
0020    if (isnan(r2)) { return true; }
0021    return r1 < r2;
0022 }
0023
0024 template<typename Real>
0025 std::array<Real, 4> polish_and_sort(Real a, Real b, Real c, Real d, Real e, std::array<Real, 4>& roots) {
0026     // Polish the roots with a Halley iterate.
0027     using std::fma;
0028     using std::abs;
0029     for (auto &r : roots) {
0030         Real df = fma(4*a, r, 3*b);
0031         df = fma(df, r, 2*c);
0032         df = fma(df, r, d);
0033         Real d2f = fma(12*a, r, 6*b);
0034         d2f = fma(d2f, r, 2*c);
0035         Real f = fma(a, r, b);
0036         f = fma(f,r,c);
0037         f = fma(f,r,d);
0038         f = fma(f,r,e);
0039         Real denom = 2*df*df - f*d2f;
0040         if (abs(denom) > (std::numeric_limits<Real>::min)())
0041         {
0042             r -= 2*f*df/denom;
0043         }
0044     }
0045     std::sort(roots.begin(), roots.end(), detail::comparator<Real>);
0046     return roots;
0047 }
0048
0049 }
0050 // Solves ax^4 + bx^3 + cx^2 + dx + e = 0.
0051 // Only returns the real roots, as these are the only roots of interest in ray intersection problems.
0052 // Follows Graphics Gems V: https://github.com/erich666/GraphicsGems/blob/master/gems/Roots3And4.c
0053 template<typename Real>
0054 std::array<Real, 4> quartic_roots(Real a, Real b, Real c, Real d, Real e) {
0055     using std::abs;
0056     using std::sqrt;
0057     auto nan = std::numeric_limits<Real>::quiet_NaN();
0058     std::array<Real, 4> roots{nan, nan, nan, nan};
0059     if (abs(a) <= (std::numeric_limits<Real>::min)()) {
0060         auto cbrts = cubic_roots(b, c, d, e);
0061         roots[0] = cbrts[0];
0062         roots[1] = cbrts[1];
0063         roots[2] = cbrts[2];
0064         if (b == 0 && c == 0 && d == 0 && e == 0) {
0065            roots[3] = 0;
0066         }
0067         return detail::polish_and_sort(a, b, c, d, e, roots);
0068     }
0069     if (abs(e) <= (std::numeric_limits<Real>::min)()) {
0070         auto v = cubic_roots(a, b, c, d);
0071         roots[0] = v[0];
0072         roots[1] = v[1];
0073         roots[2] = v[2];
0074         roots[3] = 0;
0075         return detail::polish_and_sort(a, b, c, d, e, roots);
0076     }
0077     // Now solve x^4 + Ax^3 + Bx^2 + Cx + D = 0.
0078     Real A = b/a;
0079     Real B = c/a;
0080     Real C = d/a;
0081     Real D = e/a;
0082     Real Asq = A*A;
0083     // Let x = y - A/4:
0084     // Mathematica: Expand[(y - A/4)^4 + A*(y - A/4)^3 + B*(y - A/4)^2 + C*(y - A/4) + D]
0085     // We now solve the depressed quartic y^4 + py^2 + qy + r = 0.
0086     Real p = B - 3*Asq/8;
0087     Real q = C - A*B/2 + Asq*A/8;
0088     Real r = D - A*C/4 + Asq*B/16 - 3*Asq*Asq/256;
0089     if (abs(r) <= (std::numeric_limits<Real>::min)()) {
0090         auto [r1, r2, r3] = cubic_roots(Real(1), Real(0), p, q);
0091         r1 -= A/4;
0092         r2 -= A/4;
0093         r3 -= A/4;
0094         roots[0] = r1;
0095         roots[1] = r2;
0096         roots[2] = r3;
0097         roots[3] = -A/4;
0098         return detail::polish_and_sort(a, b, c, d, e, roots);
0099     }
0100     // Biquadratic case:
0101     if (abs(q) <= (std::numeric_limits<Real>::min)()) {
0102         auto [r1, r2] = quadratic_roots(Real(1), p, r);
0103         if (r1 >= 0) {
0104            Real rtr = sqrt(r1);
0105            roots[0] = rtr - A/4;
0106            roots[1] = -rtr - A/4;
0107         }
0108         if (r2 >= 0) {
0109            Real rtr = sqrt(r2);
0110            roots[2] = rtr - A/4;
0111            roots[3] = -rtr - A/4;
0112         }
0113         return detail::polish_and_sort(a, b, c, d, e, roots);
0114     }
0115
0116     // Now split the depressed quartic into two quadratics:
0117     // y^4 + py^2 + qy + r = (y^2 + sy + u)(y^2 - sy + v) = y^4 + (v+u-s^2)y^2 + s(v - u)y + uv
0118     // So p = v+u-s^2, q = s(v - u), r = uv.
0119     // Then (v+u)^2 - (v-u)^2 = 4uv = 4r = (p+s^2)^2 - q^2/s^2.
0120     // Multiply through by s^2 to get s^2(p+s^2)^2 - q^2 - 4rs^2 = 0, which is a cubic in s^2.
0121     // Then we let z = s^2, to get
0122     // z^3 + 2pz^2 + (p^2 - 4r)z - q^2 = 0.
0123     auto z_roots = cubic_roots(Real(1), 2*p, p*p - 4*r, -q*q);
0124     // z = s^2, so s = sqrt(z).
0125     // Hence we require a root > 0, and for the sake of sanity we should take the largest one:
0126     Real largest_root = std::numeric_limits<Real>::lowest();
0127     for (auto z : z_roots) {
0128         if (z > largest_root) {
0129             largest_root = z;
0130         }
0131     }
0132     // No real roots:
0133     if (largest_root <= 0) {
0134       return roots;
0135     }
0136     Real s = sqrt(largest_root);
0137     // s is nonzero, because we took care of the biquadratic case.
0138     Real v = (p + largest_root + q/s)/2;
0139     Real u = v - q/s;
0140     // Now solve y^2 + sy + u = 0:
0141     auto [root0, root1] = quadratic_roots(Real(1), s, u);
0142
0143     // Now solve y^2 - sy + v = 0:
0144     auto [root2, root3] = quadratic_roots(Real(1), -s, v);
0145     roots[0] = root0;
0146     roots[1] = root1;
0147     roots[2] = root2;
0148     roots[3] = root3;
0149
0150     for (auto& r : roots) {
0151         r -= A/4;
0152     }
0153     return detail::polish_and_sort(a, b, c, d, e, roots);
0154 }
0155
0156 }
0157 #endif