diff --git a/include/boost/math/tools/quartic_roots.hpp b/include/boost/math/tools/quartic_roots.hpp index 56e46fd03..88a4f7253 100644 --- a/include/boost/math/tools/quartic_roots.hpp +++ b/include/boost/math/tools/quartic_roots.hpp @@ -135,7 +135,7 @@ std::array quartic_roots(Real a, Real b, Real c, Real d, Real e) { } Real s = sqrt(largest_root); // s is nonzero, because we took care of the biquadratic case. - Real v = (p + s*s + q/s)/2; + Real v = (p + largest_root + q/s)/2; Real u = v - q/s; // Now solve y^2 + sy + u = 0: auto [root0, root1] = quadratic_roots(Real(1), s, u); diff --git a/test/quartic_roots_test.cpp b/test/quartic_roots_test.cpp index 8f88dbb85..173e86ee9 100644 --- a/test/quartic_roots_test.cpp +++ b/test/quartic_roots_test.cpp @@ -142,6 +142,21 @@ void issue_825() { CHECK_NAN(roots[3]); } +void issue_1055() { + double a = 1.0; + double b = -547.5045576653938; + double c = 75042.069484941996; + double d = 273.7522788326969; + double e = 0.24965766552610175; + std::array roots = boost::math::tools::quartic_roots(a, b, c, d, e); + // This is accurate to 1e-9 on every platform *except* cygwin/g++11/c++17: + CHECK_ABSOLUTE_ERROR(-0.00182420203946279, roots[0], 1e-6); + CHECK_ABSOLUTE_ERROR(-0.00182370927680797, roots[1], 1e-6); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); +} + + int main() { test_zero_coefficients(); @@ -150,5 +165,6 @@ int main() test_zero_coefficients(); #endif issue_825(); + issue_1055(); return boost::math::test::report_errors(); }