diff options
Diffstat (limited to 'boost/math/tools/roots.hpp')
| -rw-r--r-- | boost/math/tools/roots.hpp | 22 |
1 files changed, 16 insertions, 6 deletions
diff --git a/boost/math/tools/roots.hpp b/boost/math/tools/roots.hpp index 356197dcf2..2442f5c2d1 100644 --- a/boost/math/tools/roots.hpp +++ b/boost/math/tools/roots.hpp @@ -109,13 +109,13 @@ std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, c static const char* function = "boost::math::tools::bisect<%1%>"; if(min >= max) { - policies::raise_evaluation_error(function, - "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol); + return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, + "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol)); } if(fmin * fmax >= 0) { - policies::raise_evaluation_error(function, - "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol); + return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, + "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol)); } // @@ -302,7 +302,7 @@ T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_i if(0 == f0) break; - if((f1 == 0) && (f2 == 0)) + if(f1 == 0) { // Oops zero derivative!!! #ifdef BOOST_MATH_INSTRUMENT @@ -329,8 +329,18 @@ T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_i delta = denom / num; if(delta * f1 / f0 < 0) { - // probably cancellation error, try a Newton step instead: + // Oh dear, we have a problem as Newton and Halley steps + // disagree about which way we should move. Probably + // there is cancelation error in the calculation of the + // Halley step, or else the derivatives are so small + // that their values are basically trash. We will move + // in the direction indicated by a Newton step, but + // by no more than twice the current guess value, otherwise + // we can jump way out of bounds if we're not careful. + // See https://svn.boost.org/trac/boost/ticket/8314. delta = f0 / f1; + if(fabs(delta) > 2 * fabs(guess)) + delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess); } } else |