Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: 
Chris@16: #ifndef BOOST_MATH_TOOLS_MINIMA_HPP
Chris@16: #define BOOST_MATH_TOOLS_MINIMA_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <utility>
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <boost/math/tools/precision.hpp>
Chris@16: #include <boost/math/policies/policy.hpp>
Chris@16: #include <boost/cstdint.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace tools{
Chris@16: 
Chris@16: template <class F, class T>
Chris@16: std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
Chris@16:    T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
Chris@16:    T x;  // minima so far
Chris@16:    T w;  // second best point
Chris@16:    T v;  // previous value of w
Chris@16:    T u;  // most recent evaluation point
Chris@16:    T delta;  // The distance moved in the last step
Chris@16:    T delta2; // The distance moved in the step before last
Chris@16:    T fu, fv, fw, fx;  // function evaluations at u, v, w, x
Chris@16:    T mid; // midpoint of min and max
Chris@16:    T fract1, fract2;  // minimal relative movement in x
Chris@16: 
Chris@16:    static const T golden = 0.3819660f;  // golden ratio, don't need too much precision here!
Chris@16: 
Chris@16:    x = w = v = max;
Chris@16:    fw = fv = fx = f(x);
Chris@16:    delta2 = delta = 0;
Chris@16: 
Chris@16:    uintmax_t count = max_iter;
Chris@16: 
Chris@16:    do{
Chris@16:       // get midpoint
Chris@16:       mid = (min + max) / 2;
Chris@16:       // work out if we're done already:
Chris@16:       fract1 = tolerance * fabs(x) + tolerance / 4;
Chris@16:       fract2 = 2 * fract1;
Chris@16:       if(fabs(x - mid) <= (fract2 - (max - min) / 2))
Chris@16:          break;
Chris@16: 
Chris@16:       if(fabs(delta2) > fract1)
Chris@16:       {
Chris@16:          // try and construct a parabolic fit:
Chris@16:          T r = (x - w) * (fx - fv);
Chris@16:          T q = (x - v) * (fx - fw);
Chris@16:          T p = (x - v) * q - (x - w) * r;
Chris@16:          q = 2 * (q - r);
Chris@16:          if(q > 0)
Chris@16:             p = -p;
Chris@16:          q = fabs(q);
Chris@16:          T td = delta2;
Chris@16:          delta2 = delta;
Chris@16:          // determine whether a parabolic step is acceptible or not:
Chris@16:          if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))
Chris@16:          {
Chris@16:             // nope, try golden section instead
Chris@16:             delta2 = (x >= mid) ? min - x : max - x;
Chris@16:             delta = golden * delta2;
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             // whew, parabolic fit:
Chris@16:             delta = p / q;
Chris@16:             u = x + delta;
Chris@16:             if(((u - min) < fract2) || ((max- u) < fract2))
Chris@16:                delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // golden section:
Chris@16:          delta2 = (x >= mid) ? min - x : max - x;
Chris@16:          delta = golden * delta2;
Chris@16:       }
Chris@16:       // update current position:
Chris@16:       u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
Chris@16:       fu = f(u);
Chris@16:       if(fu <= fx)
Chris@16:       {
Chris@16:          // good new point is an improvement!
Chris@16:          // update brackets:
Chris@16:          if(u >= x)
Chris@16:             min = x;
Chris@16:          else
Chris@16:             max = x;
Chris@16:          // update control points:
Chris@16:          v = w;
Chris@16:          w = x;
Chris@16:          x = u;
Chris@16:          fv = fw;
Chris@16:          fw = fx;
Chris@16:          fx = fu;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // Oh dear, point u is worse than what we have already,
Chris@16:          // even so it *must* be better than one of our endpoints:
Chris@16:          if(u < x)
Chris@16:             min = u;
Chris@16:          else
Chris@16:             max = u;
Chris@16:          if((fu <= fw) || (w == x))
Chris@16:          {
Chris@16:             // however it is at least second best:
Chris@16:             v = w;
Chris@16:             w = u;
Chris@16:             fv = fw;
Chris@16:             fw = fu;
Chris@16:          }
Chris@16:          else if((fu <= fv) || (v == x) || (v == w))
Chris@16:          {
Chris@16:             // third best:
Chris@16:             v = u;
Chris@16:             fv = fu;
Chris@16:          }
Chris@16:       }
Chris@16: 
Chris@16:    }while(--count);
Chris@16: 
Chris@16:    max_iter -= count;
Chris@16: 
Chris@16:    return std::make_pair(x, fx);
Chris@16: }
Chris@16: 
Chris@16: template <class F, class T>
Chris@16: inline std::pair<T, T> brent_find_minima(F f, T min, T max, int digits)
Chris@16: {
Chris@16:    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16:    return brent_find_minima(f, min, max, digits, m);
Chris@16: }
Chris@16: 
Chris@16: }}} // namespaces
Chris@16: 
Chris@16: #endif
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: