// (C) Copyright John Maddock 2006. | |
// Use, modification and distribution are subject to the | |
// Boost Software License, Version 1.0. (See accompanying file | |
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
#ifndef BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP | |
#define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP | |
#ifdef _MSC_VER | |
#pragma once | |
#endif | |
#include <utility> | |
#include <boost/config/no_tr1/cmath.hpp> | |
#include <stdexcept> | |
#include <boost/math/tools/config.hpp> | |
#include <boost/cstdint.hpp> | |
#include <boost/assert.hpp> | |
#include <boost/throw_exception.hpp> | |
#ifdef BOOST_MSVC | |
#pragma warning(push) | |
#pragma warning(disable: 4512) | |
#endif | |
#include <boost/math/tools/tuple.hpp> | |
#ifdef BOOST_MSVC | |
#pragma warning(pop) | |
#endif | |
#include <boost/math/special_functions/sign.hpp> | |
#include <boost/math/tools/toms748_solve.hpp> | |
#include <boost/math/policies/error_handling.hpp> | |
namespace boost{ namespace math{ namespace tools{ | |
namespace detail{ | |
template <class Tuple, class T> | |
inline void unpack_0(const Tuple& t, T& val) | |
{ val = boost::math::get<0>(t); } | |
template <class F, class T> | |
void handle_zero_derivative(F f, | |
T& last_f0, | |
const T& f0, | |
T& delta, | |
T& result, | |
T& guess, | |
const T& min, | |
const T& max) | |
{ | |
if(last_f0 == 0) | |
{ | |
// this must be the first iteration, pretend that we had a | |
// previous one at either min or max: | |
if(result == min) | |
{ | |
guess = max; | |
} | |
else | |
{ | |
guess = min; | |
} | |
unpack_0(f(guess), last_f0); | |
delta = guess - result; | |
} | |
if(sign(last_f0) * sign(f0) < 0) | |
{ | |
// we've crossed over so move in opposite direction to last step: | |
if(delta < 0) | |
{ | |
delta = (result - min) / 2; | |
} | |
else | |
{ | |
delta = (result - max) / 2; | |
} | |
} | |
else | |
{ | |
// move in same direction as last step: | |
if(delta < 0) | |
{ | |
delta = (result - max) / 2; | |
} | |
else | |
{ | |
delta = (result - min) / 2; | |
} | |
} | |
} | |
} // namespace | |
template <class F, class T, class Tol, class Policy> | |
std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) | |
{ | |
T fmin = f(min); | |
T fmax = f(max); | |
if(fmin == 0) | |
return std::make_pair(min, min); | |
if(fmax == 0) | |
return std::make_pair(max, max); | |
// | |
// Error checking: | |
// | |
static const char* function = "boost::math::tools::bisect<%1%>"; | |
if(min >= max) | |
{ | |
policies::raise_evaluation_error(function, | |
"Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol); | |
} | |
if(fmin * fmax >= 0) | |
{ | |
policies::raise_evaluation_error(function, | |
"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); | |
} | |
// | |
// Three function invocations so far: | |
// | |
boost::uintmax_t count = max_iter; | |
if(count < 3) | |
count = 0; | |
else | |
count -= 3; | |
while(count && (0 == tol(min, max))) | |
{ | |
T mid = (min + max) / 2; | |
T fmid = f(mid); | |
if((mid == max) || (mid == min)) | |
break; | |
if(fmid == 0) | |
{ | |
min = max = mid; | |
break; | |
} | |
else if(sign(fmid) * sign(fmin) < 0) | |
{ | |
max = mid; | |
fmax = fmid; | |
} | |
else | |
{ | |
min = mid; | |
fmin = fmid; | |
} | |
--count; | |
} | |
max_iter -= count; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Bisection iteration, final count = " << max_iter << std::endl; | |
static boost::uintmax_t max_count = 0; | |
if(max_iter > max_count) | |
{ | |
max_count = max_iter; | |
std::cout << "Maximum iterations: " << max_iter << std::endl; | |
} | |
#endif | |
return std::make_pair(min, max); | |
} | |
template <class F, class T, class Tol> | |
inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter) | |
{ | |
return bisect(f, min, max, tol, max_iter, policies::policy<>()); | |
} | |
template <class F, class T, class Tol> | |
inline std::pair<T, T> bisect(F f, T min, T max, Tol tol) | |
{ | |
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)(); | |
return bisect(f, min, max, tol, m, policies::policy<>()); | |
} | |
template <class F, class T> | |
T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) | |
{ | |
BOOST_MATH_STD_USING | |
T f0(0), f1, last_f0(0); | |
T result = guess; | |
T factor = static_cast<T>(ldexp(1.0, 1 - digits)); | |
T delta = 1; | |
T delta1 = tools::max_value<T>(); | |
T delta2 = tools::max_value<T>(); | |
boost::uintmax_t count(max_iter); | |
do{ | |
last_f0 = f0; | |
delta2 = delta1; | |
delta1 = delta; | |
boost::math::tie(f0, f1) = f(result); | |
if(0 == f0) | |
break; | |
if(f1 == 0) | |
{ | |
// Oops zero derivative!!! | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Newton iteration, zero derivative found" << std::endl; | |
#endif | |
detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); | |
} | |
else | |
{ | |
delta = f0 / f1; | |
} | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Newton iteration, delta = " << delta << std::endl; | |
#endif | |
if(fabs(delta * 2) > fabs(delta2)) | |
{ | |
// last two steps haven't converged, try bisection: | |
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; | |
} | |
guess = result; | |
result -= delta; | |
if(result <= min) | |
{ | |
delta = 0.5F * (guess - min); | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
else if(result >= max) | |
{ | |
delta = 0.5F * (guess - max); | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
// update brackets: | |
if(delta > 0) | |
max = guess; | |
else | |
min = guess; | |
}while(--count && (fabs(result * factor) < fabs(delta))); | |
max_iter -= count; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl; | |
static boost::uintmax_t max_count = 0; | |
if(max_iter > max_count) | |
{ | |
max_count = max_iter; | |
std::cout << "Maximum iterations: " << max_iter << std::endl; | |
} | |
#endif | |
return result; | |
} | |
template <class F, class T> | |
inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) | |
{ | |
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)(); | |
return newton_raphson_iterate(f, guess, min, max, digits, m); | |
} | |
template <class F, class T> | |
T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) | |
{ | |
BOOST_MATH_STD_USING | |
T f0(0), f1, f2; | |
T result = guess; | |
T factor = static_cast<T>(ldexp(1.0, 1 - digits)); | |
T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta | |
T last_f0 = 0; | |
T delta1 = delta; | |
T delta2 = delta; | |
bool out_of_bounds_sentry = false; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, limit = " << factor << std::endl; | |
#endif | |
boost::uintmax_t count(max_iter); | |
do{ | |
last_f0 = f0; | |
delta2 = delta1; | |
delta1 = delta; | |
boost::math::tie(f0, f1, f2) = f(result); | |
BOOST_MATH_INSTRUMENT_VARIABLE(f0); | |
BOOST_MATH_INSTRUMENT_VARIABLE(f1); | |
BOOST_MATH_INSTRUMENT_VARIABLE(f2); | |
if(0 == f0) | |
break; | |
if((f1 == 0) && (f2 == 0)) | |
{ | |
// Oops zero derivative!!! | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, zero derivative found" << std::endl; | |
#endif | |
detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); | |
} | |
else | |
{ | |
if(f2 != 0) | |
{ | |
T denom = 2 * f0; | |
T num = 2 * f1 - f0 * (f2 / f1); | |
BOOST_MATH_INSTRUMENT_VARIABLE(denom); | |
BOOST_MATH_INSTRUMENT_VARIABLE(num); | |
if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>())) | |
{ | |
// possible overflow, use Newton step: | |
delta = f0 / f1; | |
} | |
else | |
delta = denom / num; | |
if(delta * f1 / f0 < 0) | |
{ | |
// probably cancellation error, try a Newton step instead: | |
delta = f0 / f1; | |
} | |
} | |
else | |
delta = f0 / f1; | |
} | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, delta = " << delta << std::endl; | |
#endif | |
T convergence = fabs(delta / delta2); | |
if((convergence > 0.8) && (convergence < 2)) | |
{ | |
// last two steps haven't converged, try bisection: | |
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; | |
if(fabs(delta) > result) | |
delta = sign(delta) * result; // protect against huge jumps! | |
// reset delta2 so that this branch will *not* be taken on the | |
// next iteration: | |
delta2 = delta * 3; | |
BOOST_MATH_INSTRUMENT_VARIABLE(delta); | |
} | |
guess = result; | |
result -= delta; | |
BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
// check for out of bounds step: | |
if(result < min) | |
{ | |
T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min); | |
if(fabs(diff) < 1) | |
diff = 1 / diff; | |
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3)) | |
{ | |
// Only a small out of bounds step, lets assume that the result | |
// is probably approximately at min: | |
delta = 0.99f * (guess - min); | |
result = guess - delta; | |
out_of_bounds_sentry = true; // only take this branch once! | |
} | |
else | |
{ | |
delta = (guess - min) / 2; | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
} | |
else if(result > max) | |
{ | |
T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max); | |
if(fabs(diff) < 1) | |
diff = 1 / diff; | |
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3)) | |
{ | |
// Only a small out of bounds step, lets assume that the result | |
// is probably approximately at min: | |
delta = 0.99f * (guess - max); | |
result = guess - delta; | |
out_of_bounds_sentry = true; // only take this branch once! | |
} | |
else | |
{ | |
delta = (guess - max) / 2; | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
} | |
// update brackets: | |
if(delta > 0) | |
max = guess; | |
else | |
min = guess; | |
}while(--count && (fabs(result * factor) < fabs(delta))); | |
max_iter -= count; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, final count = " << max_iter << std::endl; | |
#endif | |
return result; | |
} | |
template <class F, class T> | |
inline T halley_iterate(F f, T guess, T min, T max, int digits) | |
{ | |
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)(); | |
return halley_iterate(f, guess, min, max, digits, m); | |
} | |
template <class F, class T> | |
T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) | |
{ | |
BOOST_MATH_STD_USING | |
T f0(0), f1, f2, last_f0(0); | |
T result = guess; | |
T factor = static_cast<T>(ldexp(1.0, 1 - digits)); | |
T delta = 0; | |
T delta1 = tools::max_value<T>(); | |
T delta2 = tools::max_value<T>(); | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Schroeder iteration, limit = " << factor << std::endl; | |
#endif | |
boost::uintmax_t count(max_iter); | |
do{ | |
last_f0 = f0; | |
delta2 = delta1; | |
delta1 = delta; | |
boost::math::tie(f0, f1, f2) = f(result); | |
if(0 == f0) | |
break; | |
if((f1 == 0) && (f2 == 0)) | |
{ | |
// Oops zero derivative!!! | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, zero derivative found" << std::endl; | |
#endif | |
detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); | |
} | |
else | |
{ | |
T ratio = f0 / f1; | |
if(ratio / result < 0.1) | |
{ | |
delta = ratio + (f2 / (2 * f1)) * ratio * ratio; | |
// check second derivative doesn't over compensate: | |
if(delta * ratio < 0) | |
delta = ratio; | |
} | |
else | |
delta = ratio; // fall back to Newton iteration. | |
} | |
if(fabs(delta * 2) > fabs(delta2)) | |
{ | |
// last two steps haven't converged, try bisection: | |
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; | |
} | |
guess = result; | |
result -= delta; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Halley iteration, delta = " << delta << std::endl; | |
#endif | |
if(result <= min) | |
{ | |
delta = 0.5F * (guess - min); | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
else if(result >= max) | |
{ | |
delta = 0.5F * (guess - max); | |
result = guess - delta; | |
if((result == min) || (result == max)) | |
break; | |
} | |
// update brackets: | |
if(delta > 0) | |
max = guess; | |
else | |
min = guess; | |
}while(--count && (fabs(result * factor) < fabs(delta))); | |
max_iter -= count; | |
#ifdef BOOST_MATH_INSTRUMENT | |
std::cout << "Schroeder iteration, final count = " << max_iter << std::endl; | |
static boost::uintmax_t max_count = 0; | |
if(max_iter > max_count) | |
{ | |
max_count = max_iter; | |
std::cout << "Maximum iterations: " << max_iter << std::endl; | |
} | |
#endif | |
return result; | |
} | |
template <class F, class T> | |
inline T schroeder_iterate(F f, T guess, T min, T max, int digits) | |
{ | |
boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)(); | |
return schroeder_iterate(f, guess, min, max, digits, m); | |
} | |
} // namespace tools | |
} // namespace math | |
} // namespace boost | |
#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP | |