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