158 lines · plain
1// (C) Copyright Nick Thompson 2021.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5#ifndef BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP6#define BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP7#include <array>8#include <cmath>9#include <boost/math/tools/cubic_roots.hpp>10 11namespace boost::math::tools {12 13namespace detail {14 15// Make sure the nans are always at the back of the array:16template<typename Real>17bool comparator(Real r1, Real r2) {18 using std::isnan;19 if (isnan(r1)) { return false; }20 if (isnan(r2)) { return true; }21 return r1 < r2;22}23 24template<typename Real>25std::array<Real, 4> polish_and_sort(Real a, Real b, Real c, Real d, Real e, std::array<Real, 4>& roots) {26 // Polish the roots with a Halley iterate.27 using std::fma;28 using std::abs;29 for (auto &r : roots) {30 Real df = fma(4*a, r, 3*b);31 df = fma(df, r, 2*c);32 df = fma(df, r, d);33 Real d2f = fma(12*a, r, 6*b);34 d2f = fma(d2f, r, 2*c);35 Real f = fma(a, r, b);36 f = fma(f,r,c);37 f = fma(f,r,d);38 f = fma(f,r,e);39 Real denom = 2*df*df - f*d2f;40 if (abs(denom) > (std::numeric_limits<Real>::min)())41 {42 r -= 2*f*df/denom;43 }44 }45 std::sort(roots.begin(), roots.end(), detail::comparator<Real>);46 return roots;47}48 49}50// Solves ax^4 + bx^3 + cx^2 + dx + e = 0.51// Only returns the real roots, as these are the only roots of interest in ray intersection problems.52// Follows Graphics Gems V: https://github.com/erich666/GraphicsGems/blob/master/gems/Roots3And4.c53template<typename Real>54std::array<Real, 4> quartic_roots(Real a, Real b, Real c, Real d, Real e) {55 using std::abs;56 using std::sqrt;57 auto nan = std::numeric_limits<Real>::quiet_NaN();58 std::array<Real, 4> roots{nan, nan, nan, nan};59 if (abs(a) <= (std::numeric_limits<Real>::min)()) {60 auto cbrts = cubic_roots(b, c, d, e);61 roots[0] = cbrts[0];62 roots[1] = cbrts[1];63 roots[2] = cbrts[2];64 if (b == 0 && c == 0 && d == 0 && e == 0) {65 roots[3] = 0;66 }67 return detail::polish_and_sort(a, b, c, d, e, roots);68 }69 if (abs(e) <= (std::numeric_limits<Real>::min)()) {70 auto v = cubic_roots(a, b, c, d);71 roots[0] = v[0];72 roots[1] = v[1];73 roots[2] = v[2];74 roots[3] = 0;75 return detail::polish_and_sort(a, b, c, d, e, roots);76 }77 // Now solve x^4 + Ax^3 + Bx^2 + Cx + D = 0.78 Real A = b/a;79 Real B = c/a;80 Real C = d/a;81 Real D = e/a;82 Real Asq = A*A;83 // Let x = y - A/4:84 // Mathematica: Expand[(y - A/4)^4 + A*(y - A/4)^3 + B*(y - A/4)^2 + C*(y - A/4) + D]85 // We now solve the depressed quartic y^4 + py^2 + qy + r = 0.86 Real p = B - 3*Asq/8;87 Real q = C - A*B/2 + Asq*A/8;88 Real r = D - A*C/4 + Asq*B/16 - 3*Asq*Asq/256;89 if (abs(r) <= (std::numeric_limits<Real>::min)()) {90 auto [r1, r2, r3] = cubic_roots(Real(1), Real(0), p, q);91 r1 -= A/4;92 r2 -= A/4;93 r3 -= A/4;94 roots[0] = r1;95 roots[1] = r2;96 roots[2] = r3;97 roots[3] = -A/4;98 return detail::polish_and_sort(a, b, c, d, e, roots);99 }100 // Biquadratic case:101 if (abs(q) <= (std::numeric_limits<Real>::min)()) {102 auto [r1, r2] = quadratic_roots(Real(1), p, r);103 if (r1 >= 0) {104 Real rtr = sqrt(r1);105 roots[0] = rtr - A/4;106 roots[1] = -rtr - A/4;107 }108 if (r2 >= 0) {109 Real rtr = sqrt(r2);110 roots[2] = rtr - A/4;111 roots[3] = -rtr - A/4;112 }113 return detail::polish_and_sort(a, b, c, d, e, roots);114 }115 116 // Now split the depressed quartic into two quadratics:117 // 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 + uv118 // So p = v+u-s^2, q = s(v - u), r = uv.119 // Then (v+u)^2 - (v-u)^2 = 4uv = 4r = (p+s^2)^2 - q^2/s^2.120 // 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.121 // Then we let z = s^2, to get122 // z^3 + 2pz^2 + (p^2 - 4r)z - q^2 = 0.123 auto z_roots = cubic_roots(Real(1), 2*p, p*p - 4*r, -q*q);124 // z = s^2, so s = sqrt(z).125 // Hence we require a root > 0, and for the sake of sanity we should take the largest one:126 Real largest_root = std::numeric_limits<Real>::lowest();127 for (auto z : z_roots) {128 if (z > largest_root) {129 largest_root = z;130 }131 }132 // No real roots:133 if (largest_root <= 0) {134 return roots;135 }136 Real s = sqrt(largest_root);137 // s is nonzero, because we took care of the biquadratic case.138 Real v = (p + largest_root + q/s)/2;139 Real u = v - q/s;140 // Now solve y^2 + sy + u = 0:141 auto [root0, root1] = quadratic_roots(Real(1), s, u);142 143 // Now solve y^2 - sy + v = 0:144 auto [root2, root3] = quadratic_roots(Real(1), -s, v);145 roots[0] = root0;146 roots[1] = root1;147 roots[2] = root2;148 roots[3] = root3;149 150 for (auto& r : roots) {151 r -= A/4;152 }153 return detail::polish_and_sort(a, b, c, d, e, roots);154}155 156}157#endif158