ceph/src/boost/libs/math/include/boost/math/tools/minima.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
   7 #ifndef BOOST_MATH_TOOLS_MINIMA_HPP
   8 #define BOOST_MATH_TOOLS_MINIMA_HPP
   9
  10 #ifdef _MSC_VER
  11 #pragma once
  12 #endif
  13
  14 #include <utility>
  15 #include <boost/config/no_tr1/cmath.hpp>
  16 #include <boost/math/tools/precision.hpp>
  17 #include <boost/math/policies/policy.hpp>
  18 #include <boost/cstdint.hpp>
  19
  20 namespace boost{ namespace math{ namespace tools{
  21
  22 template <class F, class T>
  23 std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter)
  24    BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
  25 {
  26    BOOST_MATH_STD_USING
  27    bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
  28    T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
  29    T x;  // minima so far
  30    T w;  // second best point
  31    T v;  // previous value of w
  32    T u;  // most recent evaluation point
  33    T delta;  // The distance moved in the last step
  34    T delta2; // The distance moved in the step before last
  35    T fu, fv, fw, fx;  // function evaluations at u, v, w, x
  36    T mid; // midpoint of min and max
  37    T fract1, fract2;  // minimal relative movement in x
  38
  39    static const T golden = 0.3819660f;  // golden ratio, don't need too much precision here!
  40
  41    x = w = v = max;
  42    fw = fv = fx = f(x);
  43    delta2 = delta = 0;
  44
  45    uintmax_t count = max_iter;
  46
  47    do{
  48       // get midpoint
  49       mid = (min + max) / 2;
  50       // work out if we're done already:
  51       fract1 = tolerance * fabs(x) + tolerance / 4;
  52       fract2 = 2 * fract1;
  53       if(fabs(x - mid) <= (fract2 - (max - min) / 2))
  54          break;
  55
  56       if(fabs(delta2) > fract1)
  57       {
  58          // try and construct a parabolic fit:
  59          T r = (x - w) * (fx - fv);
  60          T q = (x - v) * (fx - fw);
  61          T p = (x - v) * q - (x - w) * r;
  62          q = 2 * (q - r);
  63          if(q > 0)
  64             p = -p;
  65          q = fabs(q);
  66          T td = delta2;
  67          delta2 = delta;
  68          // determine whether a parabolic step is acceptible or not:
  69          if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
  70          {
  71             // nope, try golden section instead
  72             delta2 = (x >= mid) ? min - x : max - x;
  73             delta = golden * delta2;
  74          }
  75          else
  76          {
  77             // whew, parabolic fit:
  78             delta = p / q;
  79             u = x + delta;
  80             if(((u - min) < fract2) || ((max- u) < fract2))
  81                delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
  82          }
  83       }
  84       else
  85       {
  86          // golden section:
  87          delta2 = (x >= mid) ? min - x : max - x;
  88          delta = golden * delta2;
  89       }
  90       // update current position:
  91       u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
  92       fu = f(u);
  93       if(fu <= fx)
  94       {
  95          // good new point is an improvement!
  96          // update brackets:
  97          if(u >= x)
  98             min = x;
  99          else
 100             max = x;
 101          // update control points:
 102          v = w;
 103          w = x;
 104          x = u;
 105          fv = fw;
 106          fw = fx;
 107          fx = fu;
 108       }
 109       else
 110       {
 111          // Oh dear, point u is worse than what we have already,
 112          // even so it *must* be better than one of our endpoints:
 113          if(u < x)
 114             min = u;
 115          else
 116             max = u;
 117          if((fu <= fw) || (w == x))
 118          {
 119             // however it is at least second best:
 120             v = w;
 121             w = u;
 122             fv = fw;
 123             fw = fu;
 124          }
 125          else if((fu <= fv) || (v == x) || (v == w))
 126          {
 127             // third best:
 128             v = u;
 129             fv = fu;
 130          }
 131       }
 132
 133    }while(--count);
 134
 135    max_iter -= count;
 136
 137    return std::make_pair(x, fx);
 138 }
 139
 140 template <class F, class T>
 141 inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits)
 142    BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 143 {
 144    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 145    return brent_find_minima(f, min, max, digits, m);
 146 }
 147
 148 }}} // namespaces
 149
 150 #endif
 151
 152
 153
 154