Improve accuracy of quartic_roots (#1056)

boostorg · Nov 27, 2023 · e0c17f6 · e0c17f6
1 parent f7f9615
commit e0c17f6
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 1 deletion.
diff --git a/include/boost/math/tools/quartic_roots.hpp b/include/boost/math/tools/quartic_roots.hpp
@@ -135,7 +135,7 @@ std::array<Real, 4> 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
@@ -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<double, 4> roots = boost::math::tools::quartic_roots<double>(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<float>();
@@ -150,5 +165,6 @@ int main()
     test_zero_coefficients<long double>();
 #endif
     issue_825();
+    issue_1055();
     return boost::math::test::report_errors();
 }