Fix quartic roots which depressed cubic only has single real root

include/boost/math/tools/quartic_roots.hpp

@@ -122,12 +122,18 @@ std::array<Real, 4> quartic_roots(Real a, Real b, Real c, Real d, Real e) {
     // z^3 + 2pz^2 + (p^2 - 4r)z - q^2 = 0.
     auto z_roots = cubic_roots(Real(1), 2*p, p*p - 4*r, -q*q);
     // z = s^2, so s = sqrt(z).
+    // Hence we require a root > 0, and for the sake of sanity we should take the largest one:
+    Real largest_root = std::numeric_limits<Real>::lowest();
+    for (auto z : z_roots) {
+        if (z > largest_root) {
+            largest_root = z;
+	}
+    }
     // No real roots:
-    if (z_roots.back() <= 0) {
+    if (largest_root <= 0) {
       return roots;
     }
-    Real s = sqrt(z_roots.back());
-
+    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 u = v - q/s;

test/quartic_roots_test.cpp

@@ -124,6 +124,23 @@ void test_zero_coefficients()
     }
 }
 
+void issue_825() {
+    using std::sqrt;
+    using std::cbrt;
+    double a = 1;
+    double b = 1;
+    double c = 1;
+    double d = 1;
+    double e = -4;
+    auto roots = boost::math::tools::quartic_roots(a, b, c, d, e);
+    // The real roots are 1 and -1.6506
+    // Wolfram alpha: Roots[x^4 + x^3 + x^2 + x == 4]
+    double expected = (-2  - cbrt(25/(3*sqrt(6.0) - 7)) + cbrt(5*(3*sqrt(6.0) - 7)))/3;
+    CHECK_ULP_CLOSE(expected, roots[0], 5);
+    CHECK_ULP_CLOSE(1, roots[1], 5);
+    CHECK_NAN(roots[2]);
+    CHECK_NAN(roots[3]);
+}
 
 int main()
 {
@@ -132,5 +149,6 @@ int main()
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
     test_zero_coefficients<long double>();
 #endif
+    issue_825();
     return boost::math::test::report_errors();
 }