ceph/src/boost/libs/math/include/boost/math/tools/toms748_solve.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_SOLVE_ROOT_HPP
   7 #define BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/tools/precision.hpp>
  14 #include <boost/math/policies/error_handling.hpp>
  15 #include <boost/math/tools/config.hpp>
  16 #include <boost/math/special_functions/sign.hpp>
  17 #include <boost/cstdint.hpp>
  18 #include <limits>
  19
  20 #ifdef BOOST_MATH_LOG_ROOT_ITERATIONS
  21 #  define BOOST_MATH_LOGGER_INCLUDE <boost/math/tools/iteration_logger.hpp>
  22 #  include BOOST_MATH_LOGGER_INCLUDE
  23 #  undef BOOST_MATH_LOGGER_INCLUDE
  24 #else
  25 #  define BOOST_MATH_LOG_COUNT(count)
  26 #endif
  27
  28 namespace boost{ namespace math{ namespace tools{
  29
  30 template <class T>
  31 class eps_tolerance
  32 {
  33 public:
  34    eps_tolerance()
  35    {
  36       eps = 4 * tools::epsilon<T>();
  37    }
  38    eps_tolerance(unsigned bits)
  39    {
  40       BOOST_MATH_STD_USING
  41       eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon<T>()));
  42    }
  43    bool operator()(const T& a, const T& b)
  44    {
  45       BOOST_MATH_STD_USING
  46       return fabs(a - b) <= (eps * (std::min)(fabs(a), fabs(b)));
  47    }
  48 private:
  49    T eps;
  50 };
  51
  52 struct equal_floor
  53 {
  54    equal_floor(){}
  55    template <class T>
  56    bool operator()(const T& a, const T& b)
  57    {
  58       BOOST_MATH_STD_USING
  59       return floor(a) == floor(b);
  60    }
  61 };
  62
  63 struct equal_ceil
  64 {
  65    equal_ceil(){}
  66    template <class T>
  67    bool operator()(const T& a, const T& b)
  68    {
  69       BOOST_MATH_STD_USING
  70       return ceil(a) == ceil(b);
  71    }
  72 };
  73
  74 struct equal_nearest_integer
  75 {
  76    equal_nearest_integer(){}
  77    template <class T>
  78    bool operator()(const T& a, const T& b)
  79    {
  80       BOOST_MATH_STD_USING
  81       return floor(a + 0.5f) == floor(b + 0.5f);
  82    }
  83 };
  84
  85 namespace detail{
  86
  87 template <class F, class T>
  88 void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd)
  89 {
  90    //
  91    // Given a point c inside the existing enclosing interval
  92    // [a, b] sets a = c if f(c) == 0, otherwise finds the new
  93    // enclosing interval: either [a, c] or [c, b] and sets
  94    // d and fd to the point that has just been removed from
  95    // the interval.  In other words d is the third best guess
  96    // to the root.
  97    //
  98    BOOST_MATH_STD_USING  // For ADL of std math functions
  99    T tol = tools::epsilon<T>() * 2;
 100    //
 101    // If the interval [a,b] is very small, or if c is too close
 102    // to one end of the interval then we need to adjust the
 103    // location of c accordingly:
 104    //
 105    if((b - a) < 2 * tol * a)
 106    {
 107       c = a + (b - a) / 2;
 108    }
 109    else if(c <= a + fabs(a) * tol)
 110    {
 111       c = a + fabs(a) * tol;
 112    }
 113    else if(c >= b - fabs(b) * tol)
 114    {
 115       c = b - fabs(b) * tol;
 116    }
 117    //
 118    // OK, lets invoke f(c):
 119    //
 120    T fc = f(c);
 121    //
 122    // if we have a zero then we have an exact solution to the root:
 123    //
 124    if(fc == 0)
 125    {
 126       a = c;
 127       fa = 0;
 128       d = 0;
 129       fd = 0;
 130       return;
 131    }
 132    //
 133    // Non-zero fc, update the interval:
 134    //
 135    if(boost::math::sign(fa) * boost::math::sign(fc) < 0)
 136    {
 137       d = b;
 138       fd = fb;
 139       b = c;
 140       fb = fc;
 141    }
 142    else
 143    {
 144       d = a;
 145       fd = fa;
 146       a = c;
 147       fa= fc;
 148    }
 149 }
 150
 151 template <class T>
 152 inline T safe_div(T num, T denom, T r)
 153 {
 154    //
 155    // return num / denom without overflow,
 156    // return r if overflow would occur.
 157    //
 158    BOOST_MATH_STD_USING  // For ADL of std math functions
 159
 160    if(fabs(denom) < 1)
 161    {
 162       if(fabs(denom * tools::max_value<T>()) <= fabs(num))
 163          return r;
 164    }
 165    return num / denom;
 166 }
 167
 168 template <class T>
 169 inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb)
 170 {
 171    //
 172    // Performs standard secant interpolation of [a,b] given
 173    // function evaluations f(a) and f(b).  Performs a bisection
 174    // if secant interpolation would leave us very close to either
 175    // a or b.  Rationale: we only call this function when at least
 176    // one other form of interpolation has already failed, so we know
 177    // that the function is unlikely to be smooth with a root very
 178    // close to a or b.
 179    //
 180    BOOST_MATH_STD_USING  // For ADL of std math functions
 181
 182    T tol = tools::epsilon<T>() * 5;
 183    T c = a - (fa / (fb - fa)) * (b - a);
 184    if((c <= a + fabs(a) * tol) || (c >= b - fabs(b) * tol))
 185       return (a + b) / 2;
 186    return c;
 187 }
 188
 189 template <class T>
 190 T quadratic_interpolate(const T& a, const T& b, T const& d,
 191                            const T& fa, const T& fb, T const& fd,
 192                            unsigned count)
 193 {
 194    //
 195    // Performs quadratic interpolation to determine the next point,
 196    // takes count Newton steps to find the location of the
 197    // quadratic polynomial.
 198    //
 199    // Point d must lie outside of the interval [a,b], it is the third
 200    // best approximation to the root, after a and b.
 201    //
 202    // Note: this does not guarantee to find a root
 203    // inside [a, b], so we fall back to a secant step should
 204    // the result be out of range.
 205    //
 206    // Start by obtaining the coefficients of the quadratic polynomial:
 207    //
 208    T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>());
 209    T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>());
 210    A = safe_div(T(A - B), T(d - a), T(0));
 211
 212    if(A == 0)
 213    {
 214       // failure to determine coefficients, try a secant step:
 215       return secant_interpolate(a, b, fa, fb);
 216    }
 217    //
 218    // Determine the starting point of the Newton steps:
 219    //
 220    T c;
 221    if(boost::math::sign(A) * boost::math::sign(fa) > 0)
 222    {
 223       c = a;
 224    }
 225    else
 226    {
 227       c = b;
 228    }
 229    //
 230    // Take the Newton steps:
 231    //
 232    for(unsigned i = 1; i <= count; ++i)
 233    {
 234       //c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);
 235       c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a));
 236    }
 237    if((c <= a) || (c >= b))
 238    {
 239       // Oops, failure, try a secant step:
 240       c = secant_interpolate(a, b, fa, fb);
 241    }
 242    return c;
 243 }
 244
 245 template <class T>
 246 T cubic_interpolate(const T& a, const T& b, const T& d,
 247                     const T& e, const T& fa, const T& fb,
 248                     const T& fd, const T& fe)
 249 {
 250    //
 251    // Uses inverse cubic interpolation of f(x) at points
 252    // [a,b,d,e] to obtain an approximate root of f(x).
 253    // Points d and e lie outside the interval [a,b]
 254    // and are the third and forth best approximations
 255    // to the root that we have found so far.
 256    //
 257    // Note: this does not guarantee to find a root
 258    // inside [a, b], so we fall back to quadratic
 259    // interpolation in case of an erroneous result.
 260    //
 261    BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b
 262       << " d = " << d << " e = " << e << " fa = " << fa << " fb = " << fb
 263       << " fd = " << fd << " fe = " << fe);
 264    T q11 = (d - e) * fd / (fe - fd);
 265    T q21 = (b - d) * fb / (fd - fb);
 266    T q31 = (a - b) * fa / (fb - fa);
 267    T d21 = (b - d) * fd / (fd - fb);
 268    T d31 = (a - b) * fb / (fb - fa);
 269    BOOST_MATH_INSTRUMENT_CODE(
 270       "q11 = " << q11 << " q21 = " << q21 << " q31 = " << q31
 271       << " d21 = " << d21 << " d31 = " << d31);
 272    T q22 = (d21 - q11) * fb / (fe - fb);
 273    T q32 = (d31 - q21) * fa / (fd - fa);
 274    T d32 = (d31 - q21) * fd / (fd - fa);
 275    T q33 = (d32 - q22) * fa / (fe - fa);
 276    T c = q31 + q32 + q33 + a;
 277    BOOST_MATH_INSTRUMENT_CODE(
 278       "q22 = " << q22 << " q32 = " << q32 << " d32 = " << d32
 279       << " q33 = " << q33 << " c = " << c);
 280
 281    if((c <= a) || (c >= b))
 282    {
 283       // Out of bounds step, fall back to quadratic interpolation:
 284       c = quadratic_interpolate(a, b, d, fa, fb, fd, 3);
 285    BOOST_MATH_INSTRUMENT_CODE(
 286       "Out of bounds interpolation, falling back to quadratic interpolation. c = " << c);
 287    }
 288
 289    return c;
 290 }
 291
 292 } // namespace detail
 293
 294 template <class F, class T, class Tol, class Policy>
 295 std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
 296 {
 297    //
 298    // Main entry point and logic for Toms Algorithm 748
 299    // root finder.
 300    //
 301    BOOST_MATH_STD_USING  // For ADL of std math functions
 302
 303    static const char* function = "boost::math::tools::toms748_solve<%1%>";
 304
 305    boost::uintmax_t count = max_iter;
 306    T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe;
 307    static const T mu = 0.5f;
 308
 309    // initialise a, b and fa, fb:
 310    a = ax;
 311    b = bx;
 312    if(a >= b)
 313       return boost::math::detail::pair_from_single(policies::raise_domain_error(
 314          function,
 315          "Parameters a and b out of order: a=%1%", a, pol));
 316    fa = fax;
 317    fb = fbx;
 318
 319    if(tol(a, b) || (fa == 0) || (fb == 0))
 320    {
 321       max_iter = 0;
 322       if(fa == 0)
 323          b = a;
 324       else if(fb == 0)
 325          a = b;
 326       return std::make_pair(a, b);
 327    }
 328
 329    if(boost::math::sign(fa) * boost::math::sign(fb) > 0)
 330       return boost::math::detail::pair_from_single(policies::raise_domain_error(
 331          function,
 332          "Parameters a and b do not bracket the root: a=%1%", a, pol));
 333    // dummy value for fd, e and fe:
 334    fe = e = fd = 1e5F;
 335
 336    if(fa != 0)
 337    {
 338       //
 339       // On the first step we take a secant step:
 340       //
 341       c = detail::secant_interpolate(a, b, fa, fb);
 342       detail::bracket(f, a, b, c, fa, fb, d, fd);
 343       --count;
 344       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 345
 346       if(count && (fa != 0) && !tol(a, b))
 347       {
 348          //
 349          // On the second step we take a quadratic interpolation:
 350          //
 351          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
 352          e = d;
 353          fe = fd;
 354          detail::bracket(f, a, b, c, fa, fb, d, fd);
 355          --count;
 356          BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 357       }
 358    }
 359
 360    while(count && (fa != 0) && !tol(a, b))
 361    {
 362       // save our brackets:
 363       a0 = a;
 364       b0 = b;
 365       //
 366       // Starting with the third step taken
 367       // we can use either quadratic or cubic interpolation.
 368       // Cubic interpolation requires that all four function values
 369       // fa, fb, fd, and fe are distinct, should that not be the case
 370       // then variable prof will get set to true, and we'll end up
 371       // taking a quadratic step instead.
 372       //
 373       T min_diff = tools::min_value<T>() * 32;
 374       bool prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
 375       if(prof)
 376       {
 377          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
 378          BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
 379       }
 380       else
 381       {
 382          c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
 383       }
 384       //
 385       // re-bracket, and check for termination:
 386       //
 387       e = d;
 388       fe = fd;
 389       detail::bracket(f, a, b, c, fa, fb, d, fd);
 390       if((0 == --count) || (fa == 0) || tol(a, b))
 391          break;
 392       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 393       //
 394       // Now another interpolated step:
 395       //
 396       prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
 397       if(prof)
 398       {
 399          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3);
 400          BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
 401       }
 402       else
 403       {
 404          c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
 405       }
 406       //
 407       // Bracket again, and check termination condition, update e:
 408       //
 409       detail::bracket(f, a, b, c, fa, fb, d, fd);
 410       if((0 == --count) || (fa == 0) || tol(a, b))
 411          break;
 412       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 413       //
 414       // Now we take a double-length secant step:
 415       //
 416       if(fabs(fa) < fabs(fb))
 417       {
 418          u = a;
 419          fu = fa;
 420       }
 421       else
 422       {
 423          u = b;
 424          fu = fb;
 425       }
 426       c = u - 2 * (fu / (fb - fa)) * (b - a);
 427       if(fabs(c - u) > (b - a) / 2)
 428       {
 429          c = a + (b - a) / 2;
 430       }
 431       //
 432       // Bracket again, and check termination condition:
 433       //
 434       e = d;
 435       fe = fd;
 436       detail::bracket(f, a, b, c, fa, fb, d, fd);
 437       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 438       BOOST_MATH_INSTRUMENT_CODE(" tol = " << T((fabs(a) - fabs(b)) / fabs(a)));
 439       if((0 == --count) || (fa == 0) || tol(a, b))
 440          break;
 441       //
 442       // And finally... check to see if an additional bisection step is
 443       // to be taken, we do this if we're not converging fast enough:
 444       //
 445       if((b - a) < mu * (b0 - a0))
 446          continue;
 447       //
 448       // bracket again on a bisection:
 449       //
 450       e = d;
 451       fe = fd;
 452       detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd);
 453       --count;
 454       BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");
 455       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
 456    } // while loop
 457
 458    max_iter -= count;
 459    if(fa == 0)
 460    {
 461       b = a;
 462    }
 463    else if(fb == 0)
 464    {
 465       a = b;
 466    }
 467    BOOST_MATH_LOG_COUNT(max_iter)
 468    return std::make_pair(a, b);
 469 }
 470
 471 template <class F, class T, class Tol>
 472 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter)
 473 {
 474    return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>());
 475 }
 476
 477 template <class F, class T, class Tol, class Policy>
 478 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
 479 {
 480    max_iter -= 2;
 481    std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol);
 482    max_iter += 2;
 483    return r;
 484 }
 485
 486 template <class F, class T, class Tol>
 487 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter)
 488 {
 489    return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>());
 490 }
 491
 492 template <class F, class T, class Tol, class Policy>
 493 std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
 494 {
 495    BOOST_MATH_STD_USING
 496    static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>";
 497    //
 498    // Set up inital brackets:
 499    //
 500    T a = guess;
 501    T b = a;
 502    T fa = f(a);
 503    T fb = fa;
 504    //
 505    // Set up invocation count:
 506    //
 507    boost::uintmax_t count = max_iter - 1;
 508
 509    int step = 32;
 510
 511    if((fa < 0) == (guess < 0 ? !rising : rising))
 512    {
 513       //
 514       // Zero is to the right of b, so walk upwards
 515       // until we find it:
 516       //
 517       while((boost::math::sign)(fb) == (boost::math::sign)(fa))
 518       {
 519          if(count == 0)
 520             return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol));
 521          //
 522          // Heuristic: normally it's best not to increase the step sizes as we'll just end up
 523          // with a really wide range to search for the root.  However, if the initial guess was *really*
 524          // bad then we need to speed up the search otherwise we'll take forever if we're orders of
 525          // magnitude out.  This happens most often if the guess is a small value (say 1) and the result
 526          // we're looking for is close to std::numeric_limits<T>::min().
 527          //
 528          if((max_iter - count) % step == 0)
 529          {
 530             factor *= 2;
 531             if(step > 1) step /= 2;
 532          }
 533          //
 534          // Now go ahead and move our guess by "factor":
 535          //
 536          a = b;
 537          fa = fb;
 538          b *= factor;
 539          fb = f(b);
 540          --count;
 541          BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 542       }
 543    }
 544    else
 545    {
 546       //
 547       // Zero is to the left of a, so walk downwards
 548       // until we find it:
 549       //
 550       while((boost::math::sign)(fb) == (boost::math::sign)(fa))
 551       {
 552          if(fabs(a) < tools::min_value<T>())
 553          {
 554             // Escape route just in case the answer is zero!
 555             max_iter -= count;
 556             max_iter += 1;
 557             return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
 558          }
 559          if(count == 0)
 560             return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol));
 561          //
 562          // Heuristic: normally it's best not to increase the step sizes as we'll just end up
 563          // with a really wide range to search for the root.  However, if the initial guess was *really*
 564          // bad then we need to speed up the search otherwise we'll take forever if we're orders of
 565          // magnitude out.  This happens most often if the guess is a small value (say 1) and the result
 566          // we're looking for is close to std::numeric_limits<T>::min().
 567          //
 568          if((max_iter - count) % step == 0)
 569          {
 570             factor *= 2;
 571             if(step > 1) step /= 2;
 572          }
 573          //
 574          // Now go ahead and move are guess by "factor":
 575          //
 576          b = a;
 577          fb = fa;
 578          a /= factor;
 579          fa = f(a);
 580          --count;
 581          BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
 582       }
 583    }
 584    max_iter -= count;
 585    max_iter += 1;
 586    std::pair<T, T> r = toms748_solve(
 587       f,
 588       (a < 0 ? b : a),
 589       (a < 0 ? a : b),
 590       (a < 0 ? fb : fa),
 591       (a < 0 ? fa : fb),
 592       tol,
 593       count,
 594       pol);
 595    max_iter += count;
 596    BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
 597    BOOST_MATH_LOG_COUNT(max_iter)
 598    return r;
 599 }
 600
 601 template <class F, class T, class Tol>
 602 inline std::pair<T, T> bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, boost::uintmax_t& max_iter)
 603 {
 604    return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>());
 605 }
 606
 607 } // namespace tools
 608 } // namespace math
 609 } // namespace boost
 610
 611
 612 #endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
 613