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

Do not use an unguarded Newton iterate to polish roots; it goes crazy near a double root. (#759)

Browse Source
This commit is contained in:
Nick
2022-02-13 16:06:02 -08:00
committed by GitHub
parent 01a938cb11
commit d063799531
3 changed files with 125 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/roots/cubic_roots.qbk
View File

@@ -21,6 +21,10 @@ std::array<Real,3> cubic_roots(Real a, Real b, Real c, Real d);
// Recall that for a numerically computed root r satisfying r = r⁎(1+ε) for the exact root r⁎ of a function p, |p(r)| ≲ ε|rṗ(r)|.
template<typename Real>
std::array<Real, 2> cubic_root_residual(Real a, Real b, Real c, Real d, Real root);
// Computes the condition number of rootfinding. Computed via Corless, A Graduate Introduction to Numerical Methods, Section 3.2.1:
template<typename Real>
Real cubic_root_condition_number(Real a, Real b, Real c, Real d, Real root);
}
```

150
include/boost/math/tools/cubic_roots.hpp
View File

@@ -4,27 +4,27 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP
#define BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP
#include <array>
#include <algorithm>
#include <array>
#include <boost/math/special_functions/sign.hpp>
#include <boost/math/tools/roots.hpp>
namespace boost::math::tools {
// 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:
// Algorithm 954: An Accurate and Efficient Cubic and Quartic Equation Solver for Physical Applications
// However, I don't have access to that paper!
template<typename Real>
// 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: Algorithm 954: An Accurate and Efficient Cubic
// and Quartic Equation Solver for Physical Applications However, I don't have
// access to that paper!
template <typename Real>
std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
using std::sqrt;
using std::acos;
using std::cos;
using std::cbrt;
using std::abs;
using std::acos;
using std::cbrt;
using std::cos;
using std::fma;
using std::sqrt;
std::array<Real, 3> roots = {std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN()};
@@ -42,7 +42,7 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
roots[2] = 0;
return roots;
}
roots[0] = -d/c;
roots[0] = -d / c;
return roots;
}
auto [x0, x1] = quadratic_roots(b, c, d);
@@ -58,79 +58,119 @@ std::array<Real, 3> cubic_roots(Real a, Real b, Real c, Real d) {
std::sort(roots.begin(), roots.end());
return roots;
}
Real p = b/a;
Real q = c/a;
Real r = d/a;
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) {
Real p = b / a;
Real q = c / a;
Real r = d / a;
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) {
Real rtQ = sqrt(Q);
Real theta = acos(R/(Q*rtQ))/3;
Real theta = acos(R / (Q * rtQ)) / 3;
Real st = sin(theta);
Real ct = cos(theta);
roots[0] = -2*rtQ*ct - p/3;
roots[1] = -rtQ*(-ct + sqrt(Real(3))*st) - p/3;
roots[2] = rtQ*(ct + sqrt(Real(3))*st) - p/3;
roots[0] = -2 * rtQ * ct - p / 3;
roots[1] = -rtQ * (-ct + sqrt(Real(3)) * st) - p / 3;
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^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)^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));
// In Numerical Recipes, Chapter 5, Section 6, it is claimed that we
// only have one real root 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)^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));
Real B = 0;
if (A != 0) {
B = Q/A;
B = Q / A;
}
roots[0] = A + B - p/3;
roots[0] = A + B - p / 3;
// Yes, we're comparing floats for equality:
// Any perturbation pushes the roots into the complex plane; out of the bailiwick of this routine.
// Any perturbation pushes the roots into the complex plane; out of the
// bailiwick of this routine.
if (A == B || arg == 0) {
roots[1] = -A - p/3;
roots[2] = -A - p/3;
roots[1] = -A - p / 3;
roots[2] = -A - p / 3;
}
}
// Root polishing:
for (auto & r : roots) {
for (auto &r : roots) {
// Horner's method.
// Here I'll take John Gustaffson's opinion that the fma is a *distinct* operation from a*x +b:
// Make sure to compile these fmas into a single instruction and not a function call!
// (I'm looking at you Windows.)
// Here I'll take John Gustaffson's opinion that the fma is a *distinct*
// operation from a*x +b: Make sure to compile these fmas into a single
// instruction and not a function call! (I'm looking at you Windows.)
Real f = fma(a, r, b);
f = fma(f,r,c);
f = fma(f,r,d);
Real df = fma(3*a, r, 2*b);
f = fma(f, r, c);
f = fma(f, r, d);
Real df = fma(3 * a, r, 2 * b);
df = fma(df, r, c);
if (df != 0) {
// No standard library feature for fused-divide add!
r -= f/df;
Real d2f = fma(6 * a, r, 2 * b);
Real denom = 2 * df * df - f * d2f;
if (denom != 0) {
r -= 2 * f * df / denom;
} else {
r -= f / df;
}
}
}
std::sort(roots.begin(), roots.end());
return roots;
}
// 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;
// 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::abs;
using std::fma;
std::array<Real, 2> out;
Real residual = fma(a, root, b);
residual = fma(residual,root,c);
residual = fma(residual,root,d);
residual = fma(residual, root, c);
residual = fma(residual, root, d);
out[0] = residual;
Real expected_residual = fma(3*a, root, 2*b);
expected_residual = fma(expected_residual, root, c);
expected_residual = abs(root*expected_residual)*std::numeric_limits<Real>::epsilon();
out[1] = expected_residual;
// The expected residual is:
// eps*[4|ar^3| + 3|br^2| + 2|cr| + |d|]
// This can be demonstrated by assuming the coefficients and the root are
// perturbed according to the rounding model of floating point arithmetic,
// and then working through the inequalities.
root = abs(root);
Real expected_residual = fma(4 * abs(a), root, 3 * abs(b));
expected_residual = fma(expected_residual, root, 2 * abs(c));
expected_residual = fma(expected_residual, root, abs(d));
out[1] = expected_residual * std::numeric_limits<Real>::epsilon();
return out;
}
// Computes the condition number of rootfinding. This is defined in Corless, A
// Graduate Introduction to Numerical Methods, Section 3.2.1.
template <typename Real>
Real cubic_root_condition_number(Real a, Real b, Real c, Real d, Real root) {
using std::abs;
using std::fma;
// There are *absolute* condition numbers that can be defined when r = 0;
// but they basically reduce to the residual computed above.
if (root == static_cast<Real>(0)) {
return std::numeric_limits<Real>::infinity();
}
Real numerator = fma(abs(a), abs(root), abs(b));
numerator = fma(numerator, abs(root), abs(c));
numerator = fma(numerator, abs(root), abs(d));
Real denominator = fma(3 * a, root, 2 * b);
denominator = fma(denominator, root, c);
if (denominator == static_cast<Real>(0)) {
return std::numeric_limits<Real>::infinity();
}
denominator *= root;
return numerator / abs(denominator);
}
} // namespace boost::math::tools
#endif

28
test/cubic_roots_test.cpp
View File

@@ -15,6 +15,7 @@ using boost::multiprecision::float128;
using boost::math::tools::cubic_roots;
using boost::math::tools::cubic_root_residual;
using boost::math::tools::cubic_root_condition_number;
using std::cbrt;
template<class Real>
@@ -107,11 +108,34 @@ void test_zero_coefficients()
CHECK_ULP_CLOSE(r[2], roots[2], 25);
for (auto root : roots) {
auto res = cubic_root_residual(a, b,c,d, root);
CHECK_LE(abs(res[0]), 20*res[1]);
CHECK_LE(abs(res[0]), res[1]);
}
}
}
void test_ill_conditioned()
{
// An ill-conditioned root reported by SATovstun:
// "Exact" roots produced with a high-precision calcuation on Wolfram Alpha:
// NSolve[x^3 + 10000*x^2 + 200*x +1==0,x]
std::array<double, 3> expected_roots{-9999.97999997, -0.010010015026300100757327057, -0.009990014973799899662674923};
auto roots = cubic_roots<double>(1, 10000, 200, 1);
CHECK_ABSOLUTE_ERROR(expected_roots[0], roots[0], std::numeric_limits<double>::epsilon());
CHECK_ABSOLUTE_ERROR(expected_roots[1], roots[1], 1.01e-5);
CHECK_ABSOLUTE_ERROR(expected_roots[2], roots[2], 1.01e-5);
double cond = cubic_root_condition_number<double>(1, 10000, 200, 1, roots[1]);
double r1 = expected_roots[1];
// The factor of 10 is a fudge factor to make the test pass.
// Nonetheless, it does show this is basically correct:
CHECK_LE(abs(r1 - roots[1])/abs(r1), 10*std::numeric_limits<double>::epsilon()*cond);
cond = cubic_root_condition_number<double>(1, 10000, 200, 1, roots[2]);
double r2 = expected_roots[2];
// The factor of 10 is a fudge factor to make the test pass.
// Nonetheless, it does show this is basically correct:
CHECK_LE(abs(r2 - roots[2])/abs(r2), 10*std::numeric_limits<double>::epsilon()*cond);
return;
}
int main()
{
@@ -120,7 +144,7 @@ int main()
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_zero_coefficients<long double>();
#endif
test_ill_conditioned();
#ifdef BOOST_HAS_FLOAT128
// For some reason, the quadmath is way less accurate than the float/double/long double:
//test_zero_coefficients<float128>();