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 template <class Tuple, class T>
  40 inline void unpack_0(const Tuple& t, T& val)
  41 { val = boost::math::get<0>(t); }
  42
  43 template <class F, class T>
  44 void handle_zero_derivative(F f,
  45                             T& last_f0,
  46                             const T& f0,
  47                             T& delta,
  48                             T& result,
  49                             T& guess,
  50                             const T& min,
  51                             const T& max)
  52 {
  53    if(last_f0 == 0)
  54    {
  55       // this must be the first iteration, pretend that we had a
  56       // previous one at either min or max:
  57       if(result == min)
  58       {
  59          guess = max;
  60       }
  61       else
  62       {
  63          guess = min;
  64       }
  65       unpack_0(f(guess), last_f0);
  66       delta = guess - result;
  67    }
  68    if(sign(last_f0) * sign(f0) < 0)
  69    {
  70       // we've crossed over so move in opposite direction to last step:
  71       if(delta < 0)
  72       {
  73          delta = (result - min) / 2;
  74       }
  75       else
  76       {
  77          delta = (result - max) / 2;
  78       }
  79    }
  80    else
  81    {
  82       // move in same direction as last step:
  83       if(delta < 0)
  84       {
  85          delta = (result - max) / 2;
  86       }
  87       else
  88       {
  89          delta = (result - min) / 2;
  90       }
  91    }
  92 }
  93
  94 } // namespace
  95
  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)
  98 {
  99    T fmin = f(min);
 100    T fmax = f(max);
 101    if(fmin == 0)
 102       return std::make_pair(min, min);
 103    if(fmax == 0)
 104       return std::make_pair(max, max);
 105
 106    //
 107    // Error checking:
 108    //
 109    static const char* function = "boost::math::tools::bisect<%1%>";
 110    if(min >= max)
 111    {
 112       policies::raise_evaluation_error(function,
 113          "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol);
 114    }
 115    if(fmin * fmax >= 0)
 116    {
 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);
 119    }
 120
 121    //
 122    // Three function invocations so far:
 123    //
 124    boost::uintmax_t count = max_iter;
 125    if(count < 3)
 126       count = 0;
 127    else
 128       count -= 3;
 129
 130    while(count && (0 == tol(min, max)))
 131    {
 132       T mid = (min + max) / 2;
 133       T fmid = f(mid);
 134       if((mid == max) || (mid == min))
 135          break;
 136       if(fmid == 0)
 137       {
 138          min = max = mid;
 139          break;
 140       }
 141       else if(sign(fmid) * sign(fmin) < 0)
 142       {
 143          max = mid;
 144          fmax = fmid;
 145       }
 146       else
 147       {
 148          min = mid;
 149          fmin = fmid;
 150       }
 151       --count;
 152    }
 153
 154    max_iter -= count;
 155
 156 #ifdef BOOST_MATH_INSTRUMENT
 157    std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
 158
 159    static boost::uintmax_t max_count = 0;
 160    if(max_iter > max_count)
 161    {
 162       max_count = max_iter;
 163       std::cout << "Maximum iterations: " << max_iter << std::endl;
 164    }
 165 #endif
 166
 167    return std::make_pair(min, max);
 168 }
 169
 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)
 172 {
 173    return bisect(f, min, max, tol, max_iter, policies::policy<>());
 174 }
 175
 176 template <class F, class T, class Tol>
 177 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol)
 178 {
 179    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 180    return bisect(f, min, max, tol, m, policies::policy<>());
 181 }
 182
 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)
 185 {
 186    BOOST_MATH_STD_USING
 187
 188    T f0(0), f1, last_f0(0);
 189    T result = guess;
 190
 191    T factor = static_cast<T>(ldexp(1.0, 1 - digits));
 192    T delta = 1;
 193    T delta1 = tools::max_value<T>();
 194    T delta2 = tools::max_value<T>();
 195
 196    boost::uintmax_t count(max_iter);
 197
 198    do{
 199       last_f0 = f0;
 200       delta2 = delta1;
 201       delta1 = delta;
 202       boost::math::tie(f0, f1) = f(result);
 203       if(0 == f0)
 204          break;
 205       if(f1 == 0)
 206       {
 207          // Oops zero derivative!!!
 208 #ifdef BOOST_MATH_INSTRUMENT
 209          std::cout << "Newton iteration, zero derivative found" << std::endl;
 210 #endif
 211          detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 212       }
 213       else
 214       {
 215          delta = f0 / f1;
 216       }
 217 #ifdef BOOST_MATH_INSTRUMENT
 218       std::cout << "Newton iteration, delta = " << delta << std::endl;
 219 #endif
 220       if(fabs(delta * 2) > fabs(delta2))
 221       {
 222          // last two steps haven't converged, try bisection:
 223          delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
 224       }
 225       guess = result;
 226       result -= delta;
 227       if(result <= min)
 228       {
 229          delta = 0.5F * (guess - min);
 230          result = guess - delta;
 231          if((result == min) || (result == max))
 232             break;
 233       }
 234       else if(result >= max)
 235       {
 236          delta = 0.5F * (guess - max);
 237          result = guess - delta;
 238          if((result == min) || (result == max))
 239             break;
 240       }
 241       // update brackets:
 242       if(delta > 0)
 243          max = guess;
 244       else
 245          min = guess;
 246    }while(--count && (fabs(result * factor) < fabs(delta)));
 247
 248    max_iter -= count;
 249
 250 #ifdef BOOST_MATH_INSTRUMENT
 251    std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl;
 252
 253    static boost::uintmax_t max_count = 0;
 254    if(max_iter > max_count)
 255    {
 256       max_count = max_iter;
 257       std::cout << "Maximum iterations: " << max_iter << std::endl;
 258    }
 259 #endif
 260
 261    return result;
 262 }
 263
 264 template <class F, class T>
 265 inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits)
 266 {
 267    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 268    return newton_raphson_iterate(f, guess, min, max, digits, m);
 269 }
 270
 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)
 273 {
 274    BOOST_MATH_STD_USING
 275
 276    T f0(0), f1, f2;
 277    T result = guess;
 278
 279    T factor = static_cast<T>(ldexp(1.0, 1 - digits));
 280    T delta = (std::max)(T(10000000 * guess), T(10000000));  // arbitarily large delta
 281    T last_f0 = 0;
 282    T delta1 = delta;
 283    T delta2 = delta;
 284
 285    bool out_of_bounds_sentry = false;
 286
 287 #ifdef BOOST_MATH_INSTRUMENT
 288    std::cout << "Halley iteration, limit = " << factor << std::endl;
 289 #endif
 290
 291    boost::uintmax_t count(max_iter);
 292
 293    do{
 294       last_f0 = f0;
 295       delta2 = delta1;
 296       delta1 = delta;
 297       boost::math::tie(f0, f1, f2) = f(result);
 298
 299       BOOST_MATH_INSTRUMENT_VARIABLE(f0);
 300       BOOST_MATH_INSTRUMENT_VARIABLE(f1);
 301       BOOST_MATH_INSTRUMENT_VARIABLE(f2);
 302
 303       if(0 == f0)
 304          break;
 305       if((f1 == 0) && (f2 == 0))
 306       {
 307          // Oops zero derivative!!!
 308 #ifdef BOOST_MATH_INSTRUMENT
 309          std::cout << "Halley iteration, zero derivative found" << std::endl;
 310 #endif
 311          detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 312       }
 313       else
 314       {
 315          if(f2 != 0)
 316          {
 317             T denom = 2 * f0;
 318             T num = 2 * f1 - f0 * (f2 / f1);
 319
 320             BOOST_MATH_INSTRUMENT_VARIABLE(denom);
 321             BOOST_MATH_INSTRUMENT_VARIABLE(num);
 322
 323             if((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
 324             {
 325                // possible overflow, use Newton step:
 326                delta = f0 / f1;
 327             }
 328             else
 329                delta = denom / num;
 330             if(delta * f1 / f0 < 0)
 331             {
 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.
 341                delta = f0 / f1;
 342                if(fabs(delta) > 2 * fabs(guess))
 343                   delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
 344             }
 345          }
 346          else
 347             delta = f0 / f1;
 348       }
 349 #ifdef BOOST_MATH_INSTRUMENT
 350       std::cout << "Halley iteration, delta = " << delta << std::endl;
 351 #endif
 352       T convergence = fabs(delta / delta2);
 353       if((convergence > 0.8) && (convergence < 2))
 354       {
 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
 360          // next iteration:
 361          delta2 = delta * 3;
 362          BOOST_MATH_INSTRUMENT_VARIABLE(delta);
 363       }
 364       guess = result;
 365       result -= delta;
 366       BOOST_MATH_INSTRUMENT_VARIABLE(result);
 367
 368       // check for out of bounds step:
 369       if(result < min)
 370       {
 371          T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000)  : T(result / min);
 372          if(fabs(diff) < 1)
 373             diff = 1 / diff;
 374          if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 375          {
 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!
 381          }
 382          else
 383          {
 384             delta = (guess - min) / 2;
 385             result = guess - delta;
 386             if((result == min) || (result == max))
 387                break;
 388          }
 389       }
 390       else if(result > max)
 391       {
 392          T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
 393          if(fabs(diff) < 1)
 394             diff = 1 / diff;
 395          if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 396          {
 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!
 402          }
 403          else
 404          {
 405             delta = (guess - max) / 2;
 406             result = guess - delta;
 407             if((result == min) || (result == max))
 408                break;
 409          }
 410       }
 411       // update brackets:
 412       if(delta > 0)
 413          max = guess;
 414       else
 415          min = guess;
 416    }while(--count && (fabs(result * factor) < fabs(delta)));
 417
 418    max_iter -= count;
 419
 420 #ifdef BOOST_MATH_INSTRUMENT
 421    std::cout << "Halley iteration, final count = " << max_iter << std::endl;
 422 #endif
 423
 424    return result;
 425 }
 426
 427 template <class F, class T>
 428 inline T halley_iterate(F f, T guess, T min, T max, int digits)
 429 {
 430    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 431    return halley_iterate(f, guess, min, max, digits, m);
 432 }
 433
 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)
 436 {
 437    BOOST_MATH_STD_USING
 438
 439    T f0(0), f1, f2, last_f0(0);
 440    T result = guess;
 441
 442    T factor = static_cast<T>(ldexp(1.0, 1 - digits));
 443    T delta = 0;
 444    T delta1 = tools::max_value<T>();
 445    T delta2 = tools::max_value<T>();
 446
 447 #ifdef BOOST_MATH_INSTRUMENT
 448    std::cout << "Schroeder iteration, limit = " << factor << std::endl;
 449 #endif
 450
 451    boost::uintmax_t count(max_iter);
 452
 453    do{
 454       last_f0 = f0;
 455       delta2 = delta1;
 456       delta1 = delta;
 457       boost::math::tie(f0, f1, f2) = f(result);
 458       if(0 == f0)
 459          break;
 460       if((f1 == 0) && (f2 == 0))
 461       {
 462          // Oops zero derivative!!!
 463 #ifdef BOOST_MATH_INSTRUMENT
 464          std::cout << "Halley iteration, zero derivative found" << std::endl;
 465 #endif
 466          detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 467       }
 468       else
 469       {
 470          T ratio = f0 / f1;
 471          if(ratio / result < 0.1)
 472          {
 473             delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
 474             // check second derivative doesn't over compensate:
 475             if(delta * ratio < 0)
 476                delta = ratio;
 477          }
 478          else
 479             delta = ratio;  // fall back to Newton iteration.
 480       }
 481       if(fabs(delta * 2) > fabs(delta2))
 482       {
 483          // last two steps haven't converged, try bisection:
 484          delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
 485       }
 486       guess = result;
 487       result -= delta;
 488 #ifdef BOOST_MATH_INSTRUMENT
 489       std::cout << "Halley iteration, delta = " << delta << std::endl;
 490 #endif
 491       if(result <= min)
 492       {
 493          delta = 0.5F * (guess - min);
 494          result = guess - delta;
 495          if((result == min) || (result == max))
 496             break;
 497       }
 498       else if(result >= max)
 499       {
 500          delta = 0.5F * (guess - max);
 501          result = guess - delta;
 502          if((result == min) || (result == max))
 503             break;
 504       }
 505       // update brackets:
 506       if(delta > 0)
 507          max = guess;
 508       else
 509          min = guess;
 510    }while(--count && (fabs(result * factor) < fabs(delta)));
 511
 512    max_iter -= count;
 513
 514 #ifdef BOOST_MATH_INSTRUMENT
 515    std::cout << "Schroeder iteration, final count = " << max_iter << std::endl;
 516
 517    static boost::uintmax_t max_count = 0;
 518    if(max_iter > max_count)
 519    {
 520       max_count = max_iter;
 521       std::cout << "Maximum iterations: " << max_iter << std::endl;
 522    }
 523 #endif
 524
 525    return result;
 526 }
 527
 528 template <class F, class T>
 529 inline T schroeder_iterate(F f, T guess, T min, T max, int digits)
 530 {
 531    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 532    return schroeder_iterate(f, guess, min, max, digits, m);
 533 }
 534
 535 } // namespace tools
 536 } // namespace math
 537 } // namespace boost
 538
 539 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
 540
 541
 542