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/tr1/tuple.hpp>
18 #include <boost/math/tools/config.hpp>
19 #include <boost/cstdint.hpp>
20 #include <boost/assert.hpp>
21 #include <boost/throw_exception.hpp>
25 #pragma warning(disable: 4512)
27 #include <boost/tr1/tuple.hpp>
32 #include <boost/math/special_functions/sign.hpp>
33 #include <boost/math/tools/toms748_solve.hpp>
34 #include <boost/math/policies/error_handling.hpp>
36 namespace boost{ namespace math{ namespace tools{
40 template <class Tuple, class T>
41 inline void unpack_0(const Tuple& t, T& val)
42 { val = std::tr1::get<0>(t); }
44 template <class F, class T>
45 void handle_zero_derivative(F f,
56 // this must be the first iteration, pretend that we had a
57 // previous one at either min or max:
66 unpack_0(f(guess), last_f0);
67 //last_f0 = std::tr1::get<0>(f(guess));
68 delta = guess - result;
70 if(sign(last_f0) * sign(f0) < 0)
72 // we've crossed over so move in opposite direction to last step:
75 delta = (result - min) / 2;
79 delta = (result - max) / 2;
84 // move in same direction as last step:
87 delta = (result - max) / 2;
91 delta = (result - min) / 2;
98 template <class F, class T, class Tol, class Policy>
99 std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
104 return std::make_pair(min, min);
106 return std::make_pair(max, max);
111 static const char* function = "boost::math::tools::bisect<%1%>";
114 policies::raise_evaluation_error(function,
115 "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol);
119 policies::raise_evaluation_error(function,
120 "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);
124 // Three function invocations so far:
126 boost::uintmax_t count = max_iter;
132 while(count && (0 == tol(min, max)))
134 T mid = (min + max) / 2;
136 if((mid == max) || (mid == min))
143 else if(sign(fmid) * sign(fmin) < 0)
158 #ifdef BOOST_MATH_INSTRUMENT
159 std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
161 static boost::uintmax_t max_count = 0;
162 if(max_iter > max_count)
164 max_count = max_iter;
165 std::cout << "Maximum iterations: " << max_iter << std::endl;
169 return std::make_pair(min, max);
172 template <class F, class T, class Tol>
173 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)
175 return bisect(f, min, max, tol, max_iter, policies::policy<>());
178 template <class F, class T, class Tol>
179 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol)
181 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
182 return bisect(f, min, max, tol, m, policies::policy<>());
185 template <class F, class T>
186 T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
190 T f0(0), f1, last_f0(0);
193 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
195 T delta1 = tools::max_value<T>();
196 T delta2 = tools::max_value<T>();
198 boost::uintmax_t count(max_iter);
204 std::tr1::tie(f0, f1) = f(result);
209 // Oops zero derivative!!!
210 #ifdef BOOST_MATH_INSTRUMENT
211 std::cout << "Newton iteration, zero derivative found" << std::endl;
213 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
219 #ifdef BOOST_MATH_INSTRUMENT
220 std::cout << "Newton iteration, delta = " << delta << std::endl;
222 if(fabs(delta * 2) > fabs(delta2))
224 // last two steps haven't converged, try bisection:
225 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
231 delta = 0.5F * (guess - min);
232 result = guess - delta;
233 if((result == min) || (result == max))
236 else if(result >= max)
238 delta = 0.5F * (guess - max);
239 result = guess - delta;
240 if((result == min) || (result == max))
248 }while(--count && (fabs(result * factor) < fabs(delta)));
252 #ifdef BOOST_MATH_INSTRUMENT
253 std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
255 static boost::uintmax_t max_count = 0;
256 if(max_iter > max_count)
258 max_count = max_iter;
259 std::cout << "Maximum iterations: " << max_iter << std::endl;
266 template <class F, class T>
267 inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits)
269 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
270 return newton_raphson_iterate(f, guess, min, max, digits, m);
273 template <class F, class T>
274 T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
281 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
282 T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta
287 bool out_of_bounds_sentry = false;
289 #ifdef BOOST_MATH_INSTRUMENT
290 std::cout << "Halley iteration, limit = " << factor << std::endl;
293 boost::uintmax_t count(max_iter);
299 std::tr1::tie(f0, f1, f2) = f(result);
301 BOOST_MATH_INSTRUMENT_VARIABLE(f0);
302 BOOST_MATH_INSTRUMENT_VARIABLE(f1);
303 BOOST_MATH_INSTRUMENT_VARIABLE(f2);
307 if((f1 == 0) && (f2 == 0))
309 // Oops zero derivative!!!
310 #ifdef BOOST_MATH_INSTRUMENT
311 std::cout << "Halley iteration, zero derivative found" << std::endl;
313 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
320 T num = 2 * f1 - f0 * (f2 / f1);
322 BOOST_MATH_INSTRUMENT_VARIABLE(denom);
323 BOOST_MATH_INSTRUMENT_VARIABLE(num);
325 if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
327 // possible overflow, use Newton step:
332 if(delta * f1 / f0 < 0)
334 // probably cancellation error, try a Newton step instead:
341 #ifdef BOOST_MATH_INSTRUMENT
342 std::cout << "Halley iteration, delta = " << delta << std::endl;
344 T convergence = fabs(delta / delta2);
345 if((convergence > 0.8) && (convergence < 2))
347 // last two steps haven't converged, try bisection:
348 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
349 if(fabs(delta) > result)
350 delta = sign(delta) * result; // protect against huge jumps!
351 // reset delta2 so that this branch will *not* be taken on the
354 BOOST_MATH_INSTRUMENT_VARIABLE(delta);
358 BOOST_MATH_INSTRUMENT_VARIABLE(result);
360 // check for out of bounds step:
363 T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
366 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
368 // Only a small out of bounds step, lets assume that the result
369 // is probably approximately at min:
370 delta = 0.99f * (guess - min);
371 result = guess - delta;
372 out_of_bounds_sentry = true; // only take this branch once!
376 delta = (guess - min) / 2;
377 result = guess - delta;
378 if((result == min) || (result == max))
382 else if(result > max)
384 T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
387 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
389 // Only a small out of bounds step, lets assume that the result
390 // is probably approximately at min:
391 delta = 0.99f * (guess - max);
392 result = guess - delta;
393 out_of_bounds_sentry = true; // only take this branch once!
397 delta = (guess - max) / 2;
398 result = guess - delta;
399 if((result == min) || (result == max))
408 }while(--count && (fabs(result * factor) < fabs(delta)));
412 #ifdef BOOST_MATH_INSTRUMENT
413 std::cout << "Halley iteration, final count = " << max_iter << std::endl;
419 template <class F, class T>
420 inline T halley_iterate(F f, T guess, T min, T max, int digits)
422 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
423 return halley_iterate(f, guess, min, max, digits, m);
426 template <class F, class T>
427 T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
431 T f0(0), f1, f2, last_f0(0);
434 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
436 T delta1 = tools::max_value<T>();
437 T delta2 = tools::max_value<T>();
439 #ifdef BOOST_MATH_INSTRUMENT
440 std::cout << "Schroeder iteration, limit = " << factor << std::endl;
443 boost::uintmax_t count(max_iter);
449 std::tr1::tie(f0, f1, f2) = f(result);
452 if((f1 == 0) && (f2 == 0))
454 // Oops zero derivative!!!
455 #ifdef BOOST_MATH_INSTRUMENT
456 std::cout << "Halley iteration, zero derivative found" << std::endl;
458 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
463 if(ratio / result < 0.1)
465 delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
466 // check second derivative doesn't over compensate:
467 if(delta * ratio < 0)
471 delta = ratio; // fall back to Newton iteration.
473 if(fabs(delta * 2) > fabs(delta2))
475 // last two steps haven't converged, try bisection:
476 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
480 #ifdef BOOST_MATH_INSTRUMENT
481 std::cout << "Halley iteration, delta = " << delta << std::endl;
485 delta = 0.5F * (guess - min);
486 result = guess - delta;
487 if((result == min) || (result == max))
490 else if(result >= max)
492 delta = 0.5F * (guess - max);
493 result = guess - delta;
494 if((result == min) || (result == max))
502 }while(--count && (fabs(result * factor) < fabs(delta)));
506 #ifdef BOOST_MATH_INSTRUMENT
507 std::cout << "Schroeder iteration, final count = " << max_iter << std::endl;
509 static boost::uintmax_t max_count = 0;
510 if(max_iter > max_count)
512 max_count = max_iter;
513 std::cout << "Maximum iterations: " << max_iter << std::endl;
520 template <class F, class T>
521 inline T schroeder_iterate(F f, T guess, T min, T max, int digits)
523 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
524 return schroeder_iterate(f, guess, min, max, digits, m);
531 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP