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{
39 template <class Tuple, class T>
40 inline void unpack_0(const Tuple& t, T& val)
41 { val = boost::math::get<0>(t); }
43 template <class F, class T>
44 void handle_zero_derivative(F f,
55 // this must be the first iteration, pretend that we had a
56 // previous one at either min or max:
65 unpack_0(f(guess), last_f0);
66 delta = guess - result;
68 if(sign(last_f0) * sign(f0) < 0)
70 // we've crossed over so move in opposite direction to last step:
73 delta = (result - min) / 2;
77 delta = (result - max) / 2;
82 // move in same direction as last step:
85 delta = (result - max) / 2;
89 delta = (result - min) / 2;
96 template <class F, class T, class Tol, class Policy>
97 std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
102 return std::make_pair(min, min);
104 return std::make_pair(max, max);
109 static const char* function = "boost::math::tools::bisect<%1%>";
112 policies::raise_evaluation_error(function,
113 "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol);
117 policies::raise_evaluation_error(function,
118 "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);
122 // Three function invocations so far:
124 boost::uintmax_t count = max_iter;
130 while(count && (0 == tol(min, max)))
132 T mid = (min + max) / 2;
134 if((mid == max) || (mid == min))
141 else if(sign(fmid) * sign(fmin) < 0)
156 #ifdef BOOST_MATH_INSTRUMENT
157 std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
159 static boost::uintmax_t max_count = 0;
160 if(max_iter > max_count)
162 max_count = max_iter;
163 std::cout << "Maximum iterations: " << max_iter << std::endl;
167 return std::make_pair(min, max);
170 template <class F, class T, class Tol>
171 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)
173 return bisect(f, min, max, tol, max_iter, policies::policy<>());
176 template <class F, class T, class Tol>
177 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol)
179 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
180 return bisect(f, min, max, tol, m, policies::policy<>());
183 template <class F, class T>
184 T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
188 T f0(0), f1, last_f0(0);
191 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
193 T delta1 = tools::max_value<T>();
194 T delta2 = tools::max_value<T>();
196 boost::uintmax_t count(max_iter);
202 boost::math::tie(f0, f1) = f(result);
207 // Oops zero derivative!!!
208 #ifdef BOOST_MATH_INSTRUMENT
209 std::cout << "Newton iteration, zero derivative found" << std::endl;
211 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
217 #ifdef BOOST_MATH_INSTRUMENT
218 std::cout << "Newton iteration, delta = " << delta << std::endl;
220 if(fabs(delta * 2) > fabs(delta2))
222 // last two steps haven't converged, try bisection:
223 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
229 delta = 0.5F * (guess - min);
230 result = guess - delta;
231 if((result == min) || (result == max))
234 else if(result >= max)
236 delta = 0.5F * (guess - max);
237 result = guess - delta;
238 if((result == min) || (result == max))
246 }while(--count && (fabs(result * factor) < fabs(delta)));
250 #ifdef BOOST_MATH_INSTRUMENT
251 std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
253 static boost::uintmax_t max_count = 0;
254 if(max_iter > max_count)
256 max_count = max_iter;
257 std::cout << "Maximum iterations: " << max_iter << std::endl;
264 template <class F, class T>
265 inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits)
267 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
268 return newton_raphson_iterate(f, guess, min, max, digits, m);
271 template <class F, class T>
272 T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
279 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
280 T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta
285 bool out_of_bounds_sentry = false;
287 #ifdef BOOST_MATH_INSTRUMENT
288 std::cout << "Halley iteration, limit = " << factor << std::endl;
291 boost::uintmax_t count(max_iter);
297 boost::math::tie(f0, f1, f2) = f(result);
299 BOOST_MATH_INSTRUMENT_VARIABLE(f0);
300 BOOST_MATH_INSTRUMENT_VARIABLE(f1);
301 BOOST_MATH_INSTRUMENT_VARIABLE(f2);
305 if((f1 == 0) && (f2 == 0))
307 // Oops zero derivative!!!
308 #ifdef BOOST_MATH_INSTRUMENT
309 std::cout << "Halley iteration, zero derivative found" << std::endl;
311 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
318 T num = 2 * f1 - f0 * (f2 / f1);
320 BOOST_MATH_INSTRUMENT_VARIABLE(denom);
321 BOOST_MATH_INSTRUMENT_VARIABLE(num);
323 if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
325 // possible overflow, use Newton step:
330 if(delta * f1 / f0 < 0)
332 // Oh dear, we have a problem as Newton and Halley steps
333 // disagree about which way we should move. Probably
334 // there is cancelation error in the calculation of the
335 // Halley step, or else the derivatives are so small
336 // that their values are basically trash. We will move
337 // in the direction indicated by a Newton step, but
338 // by no more than twice the current guess value, otherwise
339 // we can jump way out of bounds if we're not careful.
340 // See https://svn.boost.org/trac/boost/ticket/8314.
342 if(fabs(delta) > 2 * fabs(guess))
343 delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
349 #ifdef BOOST_MATH_INSTRUMENT
350 std::cout << "Halley iteration, delta = " << delta << std::endl;
352 T convergence = fabs(delta / delta2);
353 if((convergence > 0.8) && (convergence < 2))
355 // last two steps haven't converged, try bisection:
356 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
357 if(fabs(delta) > result)
358 delta = sign(delta) * result; // protect against huge jumps!
359 // reset delta2 so that this branch will *not* be taken on the
362 BOOST_MATH_INSTRUMENT_VARIABLE(delta);
366 BOOST_MATH_INSTRUMENT_VARIABLE(result);
368 // check for out of bounds step:
371 T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
374 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
376 // Only a small out of bounds step, lets assume that the result
377 // is probably approximately at min:
378 delta = 0.99f * (guess - min);
379 result = guess - delta;
380 out_of_bounds_sentry = true; // only take this branch once!
384 delta = (guess - min) / 2;
385 result = guess - delta;
386 if((result == min) || (result == max))
390 else if(result > max)
392 T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
395 if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
397 // Only a small out of bounds step, lets assume that the result
398 // is probably approximately at min:
399 delta = 0.99f * (guess - max);
400 result = guess - delta;
401 out_of_bounds_sentry = true; // only take this branch once!
405 delta = (guess - max) / 2;
406 result = guess - delta;
407 if((result == min) || (result == max))
416 }while(--count && (fabs(result * factor) < fabs(delta)));
420 #ifdef BOOST_MATH_INSTRUMENT
421 std::cout << "Halley iteration, final count = " << max_iter << std::endl;
427 template <class F, class T>
428 inline T halley_iterate(F f, T guess, T min, T max, int digits)
430 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
431 return halley_iterate(f, guess, min, max, digits, m);
434 template <class F, class T>
435 T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter)
439 T f0(0), f1, f2, last_f0(0);
442 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
444 T delta1 = tools::max_value<T>();
445 T delta2 = tools::max_value<T>();
447 #ifdef BOOST_MATH_INSTRUMENT
448 std::cout << "Schroeder iteration, limit = " << factor << std::endl;
451 boost::uintmax_t count(max_iter);
457 boost::math::tie(f0, f1, f2) = f(result);
460 if((f1 == 0) && (f2 == 0))
462 // Oops zero derivative!!!
463 #ifdef BOOST_MATH_INSTRUMENT
464 std::cout << "Halley iteration, zero derivative found" << std::endl;
466 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
471 if(ratio / result < 0.1)
473 delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
474 // check second derivative doesn't over compensate:
475 if(delta * ratio < 0)
479 delta = ratio; // fall back to Newton iteration.
481 if(fabs(delta * 2) > fabs(delta2))
483 // last two steps haven't converged, try bisection:
484 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
488 #ifdef BOOST_MATH_INSTRUMENT
489 std::cout << "Halley iteration, delta = " << delta << std::endl;
493 delta = 0.5F * (guess - min);
494 result = guess - delta;
495 if((result == min) || (result == max))
498 else if(result >= max)
500 delta = 0.5F * (guess - max);
501 result = guess - delta;
502 if((result == min) || (result == max))
510 }while(--count && (fabs(result * factor) < fabs(delta)));
514 #ifdef BOOST_MATH_INSTRUMENT
515 std::cout << "Schroeder iteration, final count = " << max_iter << std::endl;
517 static boost::uintmax_t max_count = 0;
518 if(max_iter > max_count)
520 max_count = max_iter;
521 std::cout << "Maximum iterations: " << max_iter << std::endl;
528 template <class F, class T>
529 inline T schroeder_iterate(F f, T guess, T min, T max, int digits)
531 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
532 return schroeder_iterate(f, guess, min, max, digits, m);
539 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP