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_NEWTON_SOLVER_HPP
7 #define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
14 #include <boost/config/no_tr1/cmath.hpp>
17 #include <boost/math/tools/config.hpp>
18 #include <boost/cstdint.hpp>
19 #include <boost/assert.hpp>
20 #include <boost/throw_exception.hpp>
24 #pragma warning(disable: 4512)
26 #include <boost/math/tools/tuple.hpp>
31 #include <boost/math/special_functions/sign.hpp>
32 #include <boost/math/tools/toms748_solve.hpp>
33 #include <boost/math/policies/error_handling.hpp>
35 namespace boost{ namespace math{ namespace tools{
41 template<int n, class T>
42 typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);
45 template <class Tuple, class T>
46 void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)
49 // Use ADL to find the right overload for get:
53 template <class Tuple, class T>
54 void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)
57 // Use ADL to find the right overload for get:
63 template <class Tuple, class T>
64 inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
67 // Rely on ADL to find the correct overload of get:
71 template <class T, class U, class V>
72 inline void unpack_tuple(const std::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T)
77 template <class T, class U, class V>
78 inline void unpack_0(const std::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T)
83 template <class F, class T>
84 void handle_zero_derivative(F f,
91 const T& max) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
95 // this must be the first iteration, pretend that we had a
96 // previous one at either min or max:
105 unpack_0(f(guess), last_f0);
106 delta = guess - result;
108 if(sign(last_f0) * sign(f0) < 0)
110 // we've crossed over so move in opposite direction to last step:
113 delta = (result - min) / 2;
117 delta = (result - max) / 2;
122 // move in same direction as last step:
125 delta = (result - max) / 2;
129 delta = (result - min) / 2;
136 template <class F, class T, class Tol, class Policy>
137 std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<Policy>::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
144 return std::make_pair(min, min);
149 return std::make_pair(max, max);
155 static const char* function = "boost::math::tools::bisect<%1%>";
158 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
159 "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol));
163 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
164 "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));
168 // Three function invocations so far:
170 boost::uintmax_t count = max_iter;
176 while(count && (0 == tol(min, max)))
178 T mid = (min + max) / 2;
180 if((mid == max) || (mid == min))
187 else if(sign(fmid) * sign(fmin) < 0)
202 #ifdef BOOST_MATH_INSTRUMENT
203 std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
205 static boost::uintmax_t max_count = 0;
206 if(max_iter > max_count)
208 max_count = max_iter;
209 std::cout << "Maximum iterations: " << max_iter << std::endl;
213 return std::make_pair(min, max);
216 template <class F, class T, class Tol>
217 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
219 return bisect(f, min, max, tol, max_iter, policies::policy<>());
222 template <class F, class T, class Tol>
223 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
225 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
226 return bisect(f, min, max, tol, m, policies::policy<>());
230 template <class F, class T>
231 T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
235 T f0(0), f1, last_f0(0);
238 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
239 T delta = tools::max_value<T>();
240 T delta1 = tools::max_value<T>();
241 T delta2 = tools::max_value<T>();
243 boost::uintmax_t count(max_iter);
249 detail::unpack_tuple(f(result), f0, f1);
255 // Oops zero derivative!!!
256 #ifdef BOOST_MATH_INSTRUMENT
257 std::cout << "Newton iteration, zero derivative found" << std::endl;
259 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
265 #ifdef BOOST_MATH_INSTRUMENT
266 std::cout << "Newton iteration, delta = " << delta << std::endl;
268 if(fabs(delta * 2) > fabs(delta2))
270 // last two steps haven't converged, try bisection:
271 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
277 delta = 0.5F * (guess - min);
278 result = guess - delta;
279 if((result == min) || (result == max))
282 else if(result >= max)
284 delta = 0.5F * (guess - max);
285 result = guess - delta;
286 if((result == min) || (result == max))
294 }while(count && (fabs(result * factor) < fabs(delta)));
298 #ifdef BOOST_MATH_INSTRUMENT
299 std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
301 static boost::uintmax_t max_count = 0;
302 if(max_iter > max_count)
304 max_count = max_iter;
305 std::cout << "Maximum iterations: " << max_iter << std::endl;
312 template <class F, class T>
313 inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
315 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
316 return newton_raphson_iterate(f, guess, min, max, digits, m);
324 static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
328 T num = 2 * f1 - f0 * (f2 / f1);
331 BOOST_MATH_INSTRUMENT_VARIABLE(denom);
332 BOOST_MATH_INSTRUMENT_VARIABLE(num);
334 if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
336 // possible overflow, use Newton step:
345 template <class Stepper, class F, class T>
346 T second_order_root_finder(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
353 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
354 T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta
359 bool out_of_bounds_sentry = false;
361 #ifdef BOOST_MATH_INSTRUMENT
362 std::cout << "Second order root iteration, limit = " << factor << std::endl;
365 boost::uintmax_t count(max_iter);
371 detail::unpack_tuple(f(result), f0, f1, f2);
374 BOOST_MATH_INSTRUMENT_VARIABLE(f0);
375 BOOST_MATH_INSTRUMENT_VARIABLE(f1);
376 BOOST_MATH_INSTRUMENT_VARIABLE(f2);
382 // Oops zero derivative!!!
383 #ifdef BOOST_MATH_INSTRUMENT
384 std::cout << "Second order root iteration, zero derivative found" << std::endl;
386 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
392 delta = Stepper::step(result, f0, f1, f2);
393 if(delta * f1 / f0 < 0)
395 // Oh dear, we have a problem as Newton and Halley steps
396 // disagree about which way we should move. Probably
397 // there is cancelation error in the calculation of the
398 // Halley step, or else the derivatives are so small
399 // that their values are basically trash. We will move
400 // in the direction indicated by a Newton step, but
401 // by no more than twice the current guess value, otherwise
402 // we can jump way out of bounds if we're not careful.
403 // See https://svn.boost.org/trac/boost/ticket/8314.
405 if(fabs(delta) > 2 * fabs(guess))
406 delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
412 #ifdef BOOST_MATH_INSTRUMENT
413 std::cout << "Second order root iteration, delta = " << delta << std::endl;
415 T convergence = fabs(delta / delta2);
416 if((convergence > 0.8) && (convergence < 2))
418 // last two steps haven't converged, try bisection:
419 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
420 if(fabs(delta) > result)
421 delta = sign(delta) * result; // protect against huge jumps!
422 // reset delta2 so that this branch will *not* be taken on the
425 BOOST_MATH_INSTRUMENT_VARIABLE(delta);
429 BOOST_MATH_INSTRUMENT_VARIABLE(result);
431 // check for out of bounds step:
434 T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
437 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
439 // Only a small out of bounds step, lets assume that the result
440 // is probably approximately at min:
441 delta = 0.99f * (guess - min);
442 result = guess - delta;
443 out_of_bounds_sentry = true; // only take this branch once!
447 delta = (guess - min) / 2;
448 result = guess - delta;
449 if((result == min) || (result == max))
453 else if(result > max)
455 T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
458 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
460 // Only a small out of bounds step, lets assume that the result
461 // is probably approximately at min:
462 delta = 0.99f * (guess - max);
463 result = guess - delta;
464 out_of_bounds_sentry = true; // only take this branch once!
468 delta = (guess - max) / 2;
469 result = guess - delta;
470 if((result == min) || (result == max))
479 } while(count && (fabs(result * factor) < fabs(delta)));
483 #ifdef BOOST_MATH_INSTRUMENT
484 std::cout << "Second order root iteration, final count = " << max_iter << std::endl;
492 template <class F, class T>
493 T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
495 return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);
498 template <class F, class T>
499 inline T halley_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
501 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
502 return halley_iterate(f, guess, min, max, digits, m);
507 struct schroder_stepper
510 static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
516 delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
517 // check second derivative doesn't over compensate:
518 if(delta * ratio < 0)
522 delta = ratio; // fall back to Newton iteration.
529 template <class F, class T>
530 T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
532 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
535 template <class F, class T>
536 inline T schroder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
538 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
539 return schroder_iterate(f, guess, min, max, digits, m);
542 // These two are the old spelling of this function, retained for backwards compatibity just in case:
544 template <class F, class T>
545 T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
547 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
550 template <class F, class T>
551 inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
553 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
554 return schroder_iterate(f, guess, min, max, digits, m);
562 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP