diff options
Diffstat (limited to 'boost/math/tools/toms748_solve.hpp')
| -rw-r--r-- | boost/math/tools/toms748_solve.hpp | 50 |
1 files changed, 37 insertions, 13 deletions
diff --git a/boost/math/tools/toms748_solve.hpp b/boost/math/tools/toms748_solve.hpp index f0c42324e8..48737a821a 100644 --- a/boost/math/tools/toms748_solve.hpp +++ b/boost/math/tools/toms748_solve.hpp @@ -17,6 +17,14 @@ #include <boost/cstdint.hpp> #include <limits> +#ifdef BOOST_MATH_LOG_ROOT_ITERATIONS +# define BOOST_MATH_LOGGER_INCLUDE <boost/math/tools/iteration_logger.hpp> +# include BOOST_MATH_LOGGER_INCLUDE +# undef BOOST_MATH_LOGGER_INCLUDE +#else +# define BOOST_MATH_LOG_COUNT(count) +#endif + namespace boost{ namespace math{ namespace tools{ template <class T> @@ -197,7 +205,7 @@ T quadratic_interpolate(const T& a, const T& b, T const& d, T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>()); A = safe_div(T(A - B), T(d - a), T(0)); - if(a == 0) + if(A == 0) { // failure to determine coefficients, try a secant step: return secant_interpolate(a, b, fa, fb); @@ -298,9 +306,9 @@ std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const a = ax; b = bx; if(a >= b) - policies::raise_domain_error( + return boost::math::detail::pair_from_single(policies::raise_domain_error( function, - "Parameters a and b out of order: a=%1%", a, pol); + "Parameters a and b out of order: a=%1%", a, pol)); fa = fax; fb = fbx; @@ -315,9 +323,9 @@ std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const } if(boost::math::sign(fa) * boost::math::sign(fb) > 0) - policies::raise_domain_error( + return boost::math::detail::pair_from_single(policies::raise_domain_error( function, - "Parameters a and b do not bracket the root: a=%1%", a, pol); + "Parameters a and b do not bracket the root: a=%1%", a, pol)); // dummy value for fd, e and fe: fe = e = fd = 1e5F; @@ -452,6 +460,7 @@ std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const { a = b; } + BOOST_MATH_LOG_COUNT(max_iter) return std::make_pair(a, b); } @@ -493,6 +502,8 @@ std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool risin // boost::uintmax_t count = max_iter - 1; + int step = 32; + if((fa < 0) == (guess < 0 ? !rising : rising)) { // @@ -502,13 +513,19 @@ std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool risin while((boost::math::sign)(fb) == (boost::math::sign)(fa)) { if(count == 0) - policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); + return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol)); // - // Heuristic: every 20 iterations we double the growth factor in case the - // initial guess was *really* bad ! + // Heuristic: normally it's best not to increase the step sizes as we'll just end up + // with a really wide range to search for the root. However, if the initial guess was *really* + // bad then we need to speed up the search otherwise we'll take forever if we're orders of + // magnitude out. This happens most often if the guess is a small value (say 1) and the result + // we're looking for is close to std::numeric_limits<T>::min(). // - if((max_iter - count) % 20 == 0) + if((max_iter - count) % step == 0) + { factor *= 2; + if(step > 1) step /= 2; + } // // Now go ahead and move our guess by "factor": // @@ -536,13 +553,19 @@ std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool risin return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); } if(count == 0) - policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); + return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol)); // - // Heuristic: every 20 iterations we double the growth factor in case the - // initial guess was *really* bad ! + // Heuristic: normally it's best not to increase the step sizes as we'll just end up + // with a really wide range to search for the root. However, if the initial guess was *really* + // bad then we need to speed up the search otherwise we'll take forever if we're orders of + // magnitude out. This happens most often if the guess is a small value (say 1) and the result + // we're looking for is close to std::numeric_limits<T>::min(). // - if((max_iter - count) % 20 == 0) + if((max_iter - count) % step == 0) + { factor *= 2; + if(step > 1) step /= 2; + } // // Now go ahead and move are guess by "factor": // @@ -567,6 +590,7 @@ std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool risin pol); max_iter += count; BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); + BOOST_MATH_LOG_COUNT(max_iter) return r; } |