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