Skip to content

Commit

Permalink
Merge pull request #716 from boostorg/cubic_roots_tlc
Browse files Browse the repository at this point in the history
Remove unicode from comments and loosen up error tolerance.
  • Loading branch information
jzmaddock authored Nov 3, 2021
2 parents 535fcc2 + 15c680c commit 7108ccc
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 9 deletions.
13 changes: 6 additions & 7 deletions include/boost/math/tools/cubic_roots.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@

namespace boost::math::tools {

// Solves ax³ + bx² + cx + d = 0.
// Solves ax^3 + bx^2 + cx + d = 0.
// Only returns the real roots, as types get weird for real coefficients and complex roots.
// Follows Numerical Recipes, Chapter 5, section 6.
// NB: A better algorithm apparently exists:
Expand All @@ -29,7 +29,7 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN()};
if (a == 0) {
// bx² + cx + d = 0:
// bx^2 + cx + d = 0:
if (b == 0) {
// cx + d = 0:
if (c == 0) {
Expand Down Expand Up @@ -64,7 +64,6 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
Real Q = (p*p - 3*q)/9;
Real R = (2*p*p*p - 9*p*q + 27*r)/54;
if (R*R < Q*Q*Q) {
//std::cout << "In the R² < Q³ branch.\n";
Real rtQ = sqrt(Q);
Real theta = acos(R/(Q*rtQ))/3;
Real st = sin(theta);
Expand All @@ -74,10 +73,10 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
roots[2] = rtQ*(ct + sqrt(Real(3))*st) - p/3;
} else {
// In Numerical Recipes, Chapter 5, Section 6, it is claimed that we only have one real root
// if R² >= Q³. But this isn't true; we can even see this from equation 5.6.18.
// if R^2 >= Q^3. But this isn't true; we can even see this from equation 5.6.18.
// The condition for having three real roots is that A = B.
// It *is* the case that if we're in this branch, and we have 3 real roots, two are a double root.
// Take (x+1)²(x-2) = x³ - 3x -2 as an example. This clearly has a double root at x = -1,
// Take (x+1)^2(x-2) = x^3 - 3x -2 as an example. This clearly has a double root at x = -1,
// and it gets sent into this branch.
Real arg = R*R - Q*Q*Q;
Real A = -boost::math::sign(R)*cbrt(abs(R) + sqrt(arg));
Expand Down Expand Up @@ -113,8 +112,8 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
return roots;
}

// Computes the empirical residual p(r) (first element) and expected residual ε|rṗ(r)| (second element) for a root.
// Recall that for a numerically computed root r satisfying r = r⁎(1+ε) of a function p, |p(r)| ≲ ε|rṗ(r)|.
// Computes the empirical residual p(r) (first element) and expected residual eps*|rp'(r)| (second element) for a root.
// Recall that for a numerically computed root r satisfying r = r_0(1+eps) of a function p, |p(r)| <= eps|rp'(r)|.
template<typename Real>
std::array<Real, 2> cubic_root_residual(Real a, Real b, Real c, Real d, Real root) {
using std::fma;
Expand Down
3 changes: 1 addition & 2 deletions test/cubic_roots_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,7 @@ void test_zero_coefficients()
d = -2;
// x^3 - 2 = 0
roots = cubic_roots(a,b,c,d);
// Let's go for equality here!
CHECK_EQUAL(roots[0], cbrt(Real(2)));
CHECK_ULP_CLOSE(roots[0], cbrt(Real(2)), 2);
CHECK_NAN(roots[1]);
CHECK_NAN(roots[2]);

Expand Down

0 comments on commit 7108ccc

Please sign in to comment.