boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix quartic roots when depressed cubic only has single real root (#838)

Browse Source
This commit is contained in:
Nick
2022-10-08 21:39:08 -07:00
committed by GitHub
parent 9346271a45
commit ea9c3a27e9
2 changed files with 33 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
include/boost/math/tools/quartic_roots.hpp
View File

@@ -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;

30
test/quartic_roots_test.cpp
View File

@@ -109,13 +109,13 @@ void test_zero_coefficients()
root = static_cast<Real>(dis(gen));
}
std::sort(r.begin(), r.end());
Real a = 1;
Real b = -(r[0] + r[1] + r[2] + r[3]);
Real c = r[0]*r[1] + r[0]*r[2] + r[0]*r[3] + r[1]*r[2] + r[1]*r[3] + r[2]*r[3];
Real d = -(r[0]*r[1]*r[2] + r[0]*r[1]*r[3] + r[0]*r[2]*r[3] + r[1]*r[2]*r[3]);
Real e = r[0]*r[1]*r[2]*r[3];
a = 1;
b = -(r[0] + r[1] + r[2] + r[3]);
c = r[0]*r[1] + r[0]*r[2] + r[0]*r[3] + r[1]*r[2] + r[1]*r[3] + r[2]*r[3];
d = -(r[0]*r[1]*r[2] + r[0]*r[1]*r[3] + r[0]*r[2]*r[3] + r[1]*r[2]*r[3]);
e = r[0]*r[1]*r[2]*r[3];
auto roots = quartic_roots(a, b, c, d, e);
roots = quartic_roots(a, b, c, d, e);
// I could check the condition number here, but this is fine right?
CHECK_ULP_CLOSE(r[0], roots[0], 160);
CHECK_ULP_CLOSE(r[1], roots[1], 260);
@@ -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;
std::array<double, 4> roots = boost::math::tools::quartic_roots<double>(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.0, 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();
}