ceph/src/boost/libs/math/include/boost/math/tools/roots.hpp

   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)
   5
   6 #ifndef BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
   7 #define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <utility>
  14 #include <boost/config/no_tr1/cmath.hpp>
  15 #include <stdexcept>
  16
  17 #include <boost/math/tools/config.hpp>
  18 #include <boost/cstdint.hpp>
  19 #include <boost/assert.hpp>
  20 #include <boost/throw_exception.hpp>
  21
  22 #ifdef BOOST_MSVC
  23 #pragma warning(push)
  24 #pragma warning(disable: 4512)
  25 #endif
  26 #include <boost/math/tools/tuple.hpp>
  27 #ifdef BOOST_MSVC
  28 #pragma warning(pop)
  29 #endif
  30
  31 #include <boost/math/special_functions/sign.hpp>
  32 #include <boost/math/tools/toms748_solve.hpp>
  33 #include <boost/math/policies/error_handling.hpp>
  34
  35 namespace boost{ namespace math{ namespace tools{
  36
  37 namespace detail{
  38
  39 namespace dummy{
  40
  41    template<int n, class T>
  42    typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);
  43 }
  44
  45 template <class Tuple, class T>
  46 void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)
  47 {
  48    using dummy::get;
  49    // Use ADL to find the right overload for get:
  50    a = get<0>(t);
  51    b = get<1>(t);
  52 }
  53 template <class Tuple, class T>
  54 void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)
  55 {
  56    using dummy::get;
  57    // Use ADL to find the right overload for get:
  58    a = get<0>(t);
  59    b = get<1>(t);
  60    c = get<2>(t);
  61 }
  62
  63 template <class Tuple, class T>
  64 inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
  65 {
  66    using dummy::get;
  67    // Rely on ADL to find the correct overload of get:
  68    val = get<0>(t);
  69 }
  70
  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)
  73 {
  74    a = p.first;
  75    b = p.second;
  76 }
  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)
  79 {
  80    a = p.first;
  81 }
  82
  83 template <class F, class T>
  84 void handle_zero_derivative(F f,
  85                             T& last_f0,
  86                             const T& f0,
  87                             T& delta,
  88                             T& result,
  89                             T& guess,
  90                             const T& min,
  91                             const T& max) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
  92 {
  93    if(last_f0 == 0)
  94    {
  95       // this must be the first iteration, pretend that we had a
  96       // previous one at either min or max:
  97       if(result == min)
  98       {
  99          guess = max;
 100       }
 101       else
 102       {
 103          guess = min;
 104       }
 105       unpack_0(f(guess), last_f0);
 106       delta = guess - result;
 107    }
 108    if(sign(last_f0) * sign(f0) < 0)
 109    {
 110       // we've crossed over so move in opposite direction to last step:
 111       if(delta < 0)
 112       {
 113          delta = (result - min) / 2;
 114       }
 115       else
 116       {
 117          delta = (result - max) / 2;
 118       }
 119    }
 120    else
 121    {
 122       // move in same direction as last step:
 123       if(delta < 0)
 124       {
 125          delta = (result - max) / 2;
 126       }
 127       else
 128       {
 129          delta = (result - min) / 2;
 130       }
 131    }
 132 }
 133
 134 } // namespace
 135
 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>())))
 138 {
 139    T fmin = f(min);
 140    T fmax = f(max);
 141    if(fmin == 0)
 142    {
 143       max_iter = 2;
 144       return std::make_pair(min, min);
 145    }
 146    if(fmax == 0)
 147    {
 148       max_iter = 2;
 149       return std::make_pair(max, max);
 150    }
 151
 152    //
 153    // Error checking:
 154    //
 155    static const char* function = "boost::math::tools::bisect<%1%>";
 156    if(min >= max)
 157    {
 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));
 160    }
 161    if(fmin * fmax >= 0)
 162    {
 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));
 165    }
 166
 167    //
 168    // Three function invocations so far:
 169    //
 170    boost::uintmax_t count = max_iter;
 171    if(count < 3)
 172       count = 0;
 173    else
 174       count -= 3;
 175
 176    while(count && (0 == tol(min, max)))
 177    {
 178       T mid = (min + max) / 2;
 179       T fmid = f(mid);
 180       if((mid == max) || (mid == min))
 181          break;
 182       if(fmid == 0)
 183       {
 184          min = max = mid;
 185          break;
 186       }
 187       else if(sign(fmid) * sign(fmin) < 0)
 188       {
 189          max = mid;
 190          fmax = fmid;
 191       }
 192       else
 193       {
 194          min = mid;
 195          fmin = fmid;
 196       }
 197       --count;
 198    }
 199
 200    max_iter -= count;
 201
 202 #ifdef BOOST_MATH_INSTRUMENT
 203    std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
 204
 205    static boost::uintmax_t max_count = 0;
 206    if(max_iter > max_count)
 207    {
 208       max_count = max_iter;
 209       std::cout << "Maximum iterations: " << max_iter << std::endl;
 210    }
 211 #endif
 212
 213    return std::make_pair(min, max);
 214 }
 215
 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>())))
 218 {
 219    return bisect(f, min, max, tol, max_iter, policies::policy<>());
 220 }
 221
 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>())))
 224 {
 225    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 226    return bisect(f, min, max, tol, m, policies::policy<>());
 227 }
 228
 229
 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>())))
 232 {
 233    BOOST_MATH_STD_USING
 234
 235    T f0(0), f1, last_f0(0);
 236    T result = guess;
 237
 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>();
 242
 243    boost::uintmax_t count(max_iter);
 244
 245    do{
 246       last_f0 = f0;
 247       delta2 = delta1;
 248       delta1 = delta;
 249       detail::unpack_tuple(f(result), f0, f1);
 250       --count;
 251       if(0 == f0)
 252          break;
 253       if(f1 == 0)
 254       {
 255          // Oops zero derivative!!!
 256 #ifdef BOOST_MATH_INSTRUMENT
 257          std::cout << "Newton iteration, zero derivative found" << std::endl;
 258 #endif
 259          detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 260       }
 261       else
 262       {
 263          delta = f0 / f1;
 264       }
 265 #ifdef BOOST_MATH_INSTRUMENT
 266       std::cout << "Newton iteration, delta = " << delta << std::endl;
 267 #endif
 268       if(fabs(delta * 2) > fabs(delta2))
 269       {
 270          // last two steps haven't converged, try bisection:
 271          delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
 272       }
 273       guess = result;
 274       result -= delta;
 275       if(result <= min)
 276       {
 277          delta = 0.5F * (guess - min);
 278          result = guess - delta;
 279          if((result == min) || (result == max))
 280             break;
 281       }
 282       else if(result >= max)
 283       {
 284          delta = 0.5F * (guess - max);
 285          result = guess - delta;
 286          if((result == min) || (result == max))
 287             break;
 288       }
 289       // update brackets:
 290       if(delta > 0)
 291          max = guess;
 292       else
 293          min = guess;
 294    }while(count && (fabs(result * factor) < fabs(delta)));
 295
 296    max_iter -= count;
 297
 298 #ifdef BOOST_MATH_INSTRUMENT
 299    std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
 300
 301    static boost::uintmax_t max_count = 0;
 302    if(max_iter > max_count)
 303    {
 304       max_count = max_iter;
 305       std::cout << "Maximum iterations: " << max_iter << std::endl;
 306    }
 307 #endif
 308
 309    return result;
 310 }
 311
 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>())))
 314 {
 315    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 316    return newton_raphson_iterate(f, guess, min, max, digits, m);
 317 }
 318
 319 namespace detail{
 320
 321    struct halley_step
 322    {
 323       template <class T>
 324       static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
 325       {
 326          using std::fabs;
 327          T denom = 2 * f0;
 328          T num = 2 * f1 - f0 * (f2 / f1);
 329          T delta;
 330
 331          BOOST_MATH_INSTRUMENT_VARIABLE(denom);
 332          BOOST_MATH_INSTRUMENT_VARIABLE(num);
 333
 334          if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
 335          {
 336             // possible overflow, use Newton step:
 337             delta = f0 / f1;
 338          }
 339          else
 340             delta = denom / num;
 341          return delta;
 342       }
 343    };
 344
 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>())))
 347    {
 348       BOOST_MATH_STD_USING
 349
 350          T f0(0), f1, f2;
 351       T result = guess;
 352
 353       T factor = static_cast<T>(ldexp(1.0, 1 - digits));
 354       T delta = (std::max)(T(10000000 * guess), T(10000000));  // arbitarily large delta
 355       T last_f0 = 0;
 356       T delta1 = delta;
 357       T delta2 = delta;
 358
 359       bool out_of_bounds_sentry = false;
 360
 361 #ifdef BOOST_MATH_INSTRUMENT
 362       std::cout << "Second order root iteration, limit = " << factor << std::endl;
 363 #endif
 364
 365       boost::uintmax_t count(max_iter);
 366
 367       do{
 368          last_f0 = f0;
 369          delta2 = delta1;
 370          delta1 = delta;
 371          detail::unpack_tuple(f(result), f0, f1, f2);
 372          --count;
 373
 374          BOOST_MATH_INSTRUMENT_VARIABLE(f0);
 375          BOOST_MATH_INSTRUMENT_VARIABLE(f1);
 376          BOOST_MATH_INSTRUMENT_VARIABLE(f2);
 377
 378          if(0 == f0)
 379             break;
 380          if(f1 == 0)
 381          {
 382             // Oops zero derivative!!!
 383 #ifdef BOOST_MATH_INSTRUMENT
 384             std::cout << "Second order root iteration, zero derivative found" << std::endl;
 385 #endif
 386             detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 387          }
 388          else
 389          {
 390             if(f2 != 0)
 391             {
 392                delta = Stepper::step(result, f0, f1, f2);
 393                if(delta * f1 / f0 < 0)
 394                {
 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.
 404                   delta = f0 / f1;
 405                   if(fabs(delta) > 2 * fabs(guess))
 406                      delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
 407                }
 408             }
 409             else
 410                delta = f0 / f1;
 411          }
 412 #ifdef BOOST_MATH_INSTRUMENT
 413          std::cout << "Second order root iteration, delta = " << delta << std::endl;
 414 #endif
 415          T convergence = fabs(delta / delta2);
 416          if((convergence > 0.8) && (convergence < 2))
 417          {
 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
 423             // next iteration:
 424             delta2 = delta * 3;
 425             BOOST_MATH_INSTRUMENT_VARIABLE(delta);
 426          }
 427          guess = result;
 428          result -= delta;
 429          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 430
 431          // check for out of bounds step:
 432          if(result < min)
 433          {
 434             T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
 435             if(fabs(diff) < 1)
 436                diff = 1 / diff;
 437             if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 438             {
 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!
 444             }
 445             else
 446             {
 447                delta = (guess - min) / 2;
 448                result = guess - delta;
 449                if((result == min) || (result == max))
 450                   break;
 451             }
 452          }
 453          else if(result > max)
 454          {
 455             T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
 456             if(fabs(diff) < 1)
 457                diff = 1 / diff;
 458             if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 459             {
 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!
 465             }
 466             else
 467             {
 468                delta = (guess - max) / 2;
 469                result = guess - delta;
 470                if((result == min) || (result == max))
 471                   break;
 472             }
 473          }
 474          // update brackets:
 475          if(delta > 0)
 476             max = guess;
 477          else
 478             min = guess;
 479       } while(count && (fabs(result * factor) < fabs(delta)));
 480
 481       max_iter -= count;
 482
 483 #ifdef BOOST_MATH_INSTRUMENT
 484       std::cout << "Second order root iteration, final count = " << max_iter << std::endl;
 485 #endif
 486
 487       return result;
 488    }
 489
 490 }
 491
 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>())))
 494 {
 495    return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);
 496 }
 497
 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>())))
 500 {
 501    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 502    return halley_iterate(f, guess, min, max, digits, m);
 503 }
 504
 505 namespace detail{
 506
 507    struct schroder_stepper
 508    {
 509       template <class T>
 510       static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
 511       {
 512          T ratio = f0 / f1;
 513          T delta;
 514          if(ratio / x < 0.1)
 515          {
 516             delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
 517             // check second derivative doesn't over compensate:
 518             if(delta * ratio < 0)
 519                delta = ratio;
 520          }
 521          else
 522             delta = ratio;  // fall back to Newton iteration.
 523          return delta;
 524       }
 525    };
 526
 527 }
 528
 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>())))
 531 {
 532    return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
 533 }
 534
 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>())))
 537 {
 538    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 539    return schroder_iterate(f, guess, min, max, digits, m);
 540 }
 541 //
 542 // These two are the old spelling of this function, retained for backwards compatibity just in case:
 543 //
 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>())))
 546 {
 547    return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
 548 }
 549
 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>())))
 552 {
 553    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 554    return schroder_iterate(f, guess, min, max, digits, m);
 555 }
 556
 557
 558 } // namespace tools
 559 } // namespace math
 560 } // namespace boost
 561
 562 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
 563