1 // (C) Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 #ifndef BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
7 #define BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
13 #include <boost/math/tools/precision.hpp>
14 #include <boost/math/policies/error_handling.hpp>
15 #include <boost/math/tools/config.hpp>
16 #include <boost/math/special_functions/sign.hpp>
17 #include <boost/cstdint.hpp>
20 namespace boost{ namespace math{ namespace tools{
26 eps_tolerance(unsigned bits)
29 eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(2 * tools::epsilon<T>()));
31 bool operator()(const T& a, const T& b)
34 return (fabs(a - b) / (std::min)(fabs(a), fabs(b))) <= eps;
44 bool operator()(const T& a, const T& b)
47 return floor(a) == floor(b);
55 bool operator()(const T& a, const T& b)
58 return ceil(a) == ceil(b);
62 struct equal_nearest_integer
64 equal_nearest_integer(){}
66 bool operator()(const T& a, const T& b)
69 return floor(a + 0.5f) == floor(b + 0.5f);
75 template <class F, class T>
76 void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd)
79 // Given a point c inside the existing enclosing interval
80 // [a, b] sets a = c if f(c) == 0, otherwise finds the new
81 // enclosing interval: either [a, c] or [c, b] and sets
82 // d and fd to the point that has just been removed from
83 // the interval. In other words d is the third best guess
86 BOOST_MATH_STD_USING // For ADL of std math functions
87 T tol = tools::epsilon<T>() * 2;
89 // If the interval [a,b] is very small, or if c is too close
90 // to one end of the interval then we need to adjust the
91 // location of c accordingly:
93 if((b - a) < 2 * tol * a)
97 else if(c <= a + fabs(a) * tol)
99 c = a + fabs(a) * tol;
101 else if(c >= b - fabs(b) * tol)
103 c = b - fabs(a) * tol;
106 // OK, lets invoke f(c):
110 // if we have a zero then we have an exact solution to the root:
121 // Non-zero fc, update the interval:
123 if(boost::math::sign(fa) * boost::math::sign(fc) < 0)
140 inline T safe_div(T num, T denom, T r)
143 // return num / denom without overflow,
144 // return r if overflow would occur.
146 BOOST_MATH_STD_USING // For ADL of std math functions
150 if(fabs(denom * tools::max_value<T>()) <= fabs(num))
157 inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb)
160 // Performs standard secant interpolation of [a,b] given
161 // function evaluations f(a) and f(b). Performs a bisection
162 // if secant interpolation would leave us very close to either
163 // a or b. Rationale: we only call this function when at least
164 // one other form of interpolation has already failed, so we know
165 // that the function is unlikely to be smooth with a root very
168 BOOST_MATH_STD_USING // For ADL of std math functions
170 T tol = tools::epsilon<T>() * 5;
171 T c = a - (fa / (fb - fa)) * (b - a);
172 if((c <= a + fabs(a) * tol) || (c >= b - fabs(b) * tol))
178 T quadratic_interpolate(const T& a, const T& b, T const& d,
179 const T& fa, const T& fb, T const& fd,
183 // Performs quadratic interpolation to determine the next point,
184 // takes count Newton steps to find the location of the
185 // quadratic polynomial.
187 // Point d must lie outside of the interval [a,b], it is the third
188 // best approximation to the root, after a and b.
190 // Note: this does not guarentee to find a root
191 // inside [a, b], so we fall back to a secant step should
192 // the result be out of range.
194 // Start by obtaining the coefficients of the quadratic polynomial:
196 T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>());
197 T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>());
198 A = safe_div(T(A - B), T(d - a), T(0));
202 // failure to determine coefficients, try a secant step:
203 return secant_interpolate(a, b, fa, fb);
206 // Determine the starting point of the Newton steps:
209 if(boost::math::sign(A) * boost::math::sign(fa) > 0)
218 // Take the Newton steps:
220 for(unsigned i = 1; i <= count; ++i)
222 //c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);
223 c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a));
225 if((c <= a) || (c >= b))
227 // Oops, failure, try a secant step:
228 c = secant_interpolate(a, b, fa, fb);
234 T cubic_interpolate(const T& a, const T& b, const T& d,
235 const T& e, const T& fa, const T& fb,
236 const T& fd, const T& fe)
239 // Uses inverse cubic interpolation of f(x) at points
240 // [a,b,d,e] to obtain an approximate root of f(x).
241 // Points d and e lie outside the interval [a,b]
242 // and are the third and forth best approximations
243 // to the root that we have found so far.
245 // Note: this does not guarentee to find a root
246 // inside [a, b], so we fall back to quadratic
247 // interpolation in case of an erroneous result.
249 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b
250 << " d = " << d << " e = " << e << " fa = " << fa << " fb = " << fb
251 << " fd = " << fd << " fe = " << fe);
252 T q11 = (d - e) * fd / (fe - fd);
253 T q21 = (b - d) * fb / (fd - fb);
254 T q31 = (a - b) * fa / (fb - fa);
255 T d21 = (b - d) * fd / (fd - fb);
256 T d31 = (a - b) * fb / (fb - fa);
257 BOOST_MATH_INSTRUMENT_CODE(
258 "q11 = " << q11 << " q21 = " << q21 << " q31 = " << q31
259 << " d21 = " << d21 << " d31 = " << d31);
260 T q22 = (d21 - q11) * fb / (fe - fb);
261 T q32 = (d31 - q21) * fa / (fd - fa);
262 T d32 = (d31 - q21) * fd / (fd - fa);
263 T q33 = (d32 - q22) * fa / (fe - fa);
264 T c = q31 + q32 + q33 + a;
265 BOOST_MATH_INSTRUMENT_CODE(
266 "q22 = " << q22 << " q32 = " << q32 << " d32 = " << d32
267 << " q33 = " << q33 << " c = " << c);
269 if((c <= a) || (c >= b))
271 // Out of bounds step, fall back to quadratic interpolation:
272 c = quadratic_interpolate(a, b, d, fa, fb, fd, 3);
273 BOOST_MATH_INSTRUMENT_CODE(
274 "Out of bounds interpolation, falling back to quadratic interpolation. c = " << c);
280 } // namespace detail
282 template <class F, class T, class Tol, class Policy>
283 std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
286 // Main entry point and logic for Toms Algorithm 748
289 BOOST_MATH_STD_USING // For ADL of std math functions
291 static const char* function = "boost::math::tools::toms748_solve<%1%>";
293 boost::uintmax_t count = max_iter;
294 T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe;
295 static const T mu = 0.5f;
297 // initialise a, b and fa, fb:
301 policies::raise_domain_error(
303 "Parameters a and b out of order: a=%1%", a, pol);
307 if(tol(a, b) || (fa == 0) || (fb == 0))
314 return std::make_pair(a, b);
317 if(boost::math::sign(fa) * boost::math::sign(fb) > 0)
318 policies::raise_domain_error(
320 "Parameters a and b do not bracket the root: a=%1%", a, pol);
321 // dummy value for fd, e and fe:
327 // On the first step we take a secant step:
329 c = detail::secant_interpolate(a, b, fa, fb);
330 detail::bracket(f, a, b, c, fa, fb, d, fd);
332 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
334 if(count && (fa != 0) && !tol(a, b))
337 // On the second step we take a quadratic interpolation:
339 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
342 detail::bracket(f, a, b, c, fa, fb, d, fd);
344 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
348 while(count && (fa != 0) && !tol(a, b))
350 // save our brackets:
354 // Starting with the third step taken
355 // we can use either quadratic or cubic interpolation.
356 // Cubic interpolation requires that all four function values
357 // fa, fb, fd, and fe are distinct, should that not be the case
358 // then variable prof will get set to true, and we'll end up
359 // taking a quadratic step instead.
361 T min_diff = tools::min_value<T>() * 32;
362 bool prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
365 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
366 BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
370 c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
373 // re-bracket, and check for termination:
377 detail::bracket(f, a, b, c, fa, fb, d, fd);
378 if((0 == --count) || (fa == 0) || tol(a, b))
380 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
382 // Now another interpolated step:
384 prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
387 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3);
388 BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
392 c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
395 // Bracket again, and check termination condition, update e:
397 detail::bracket(f, a, b, c, fa, fb, d, fd);
398 if((0 == --count) || (fa == 0) || tol(a, b))
400 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
402 // Now we take a double-length secant step:
404 if(fabs(fa) < fabs(fb))
414 c = u - 2 * (fu / (fb - fa)) * (b - a);
415 if(fabs(c - u) > (b - a) / 2)
420 // Bracket again, and check termination condition:
424 detail::bracket(f, a, b, c, fa, fb, d, fd);
425 if((0 == --count) || (fa == 0) || tol(a, b))
427 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
429 // And finally... check to see if an additional bisection step is
430 // to be taken, we do this if we're not converging fast enough:
432 if((b - a) < mu * (b0 - a0))
435 // bracket again on a bisection:
439 detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd);
441 BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");
442 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
454 return std::make_pair(a, b);
457 template <class F, class T, class Tol>
458 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter)
460 return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>());
463 template <class F, class T, class Tol, class Policy>
464 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
467 std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol);
472 template <class F, class T, class Tol>
473 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter)
475 return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>());
478 template <class F, class T, class Tol, class Policy>
479 std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
482 static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>";
484 // Set up inital brackets:
491 // Set up invocation count:
493 boost::uintmax_t count = max_iter - 1;
495 if((fa < 0) == (guess < 0 ? !rising : rising))
498 // Zero is to the right of b, so walk upwards
501 while((boost::math::sign)(fb) == (boost::math::sign)(fa))
504 policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol);
506 // Heuristic: every 20 iterations we double the growth factor in case the
507 // initial guess was *really* bad !
509 if((max_iter - count) % 20 == 0)
512 // Now go ahead and move our guess by "factor":
519 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
525 // Zero is to the left of a, so walk downwards
528 while((boost::math::sign)(fb) == (boost::math::sign)(fa))
530 if(fabs(a) < tools::min_value<T>())
532 // Escape route just in case the answer is zero!
535 return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
538 policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol);
540 // Heuristic: every 20 iterations we double the growth factor in case the
541 // initial guess was *really* bad !
543 if((max_iter - count) % 20 == 0)
546 // Now go ahead and move are guess by "factor":
553 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
558 std::pair<T, T> r = toms748_solve(
568 BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
572 template <class F, class T, class Tol>
573 inline std::pair<T, T> bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, boost::uintmax_t& max_iter)
575 return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>());
583 #endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP