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)
7 #ifndef BOOST_MATH_TOOLS_MINIMA_HPP
8 #define BOOST_MATH_TOOLS_MINIMA_HPP
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>
20 namespace boost{ namespace math{ namespace tools{
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>())))
27 bits = (std::min)(policies::digits<T, policies::policy<> >() / 2, bits);
28 T tolerance = static_cast<T>(ldexp(1.0, 1-bits));
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
39 static const T golden = 0.3819660f; // golden ratio, don't need too much precision here!
45 uintmax_t count = max_iter;
49 mid = (min + max) / 2;
50 // work out if we're done already:
51 fract1 = tolerance * fabs(x) + tolerance / 4;
53 if(fabs(x - mid) <= (fract2 - (max - min) / 2))
56 if(fabs(delta2) > fract1)
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;
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)))
71 // nope, try golden section instead
72 delta2 = (x >= mid) ? min - x : max - x;
73 delta = golden * delta2;
77 // whew, parabolic fit:
80 if(((u - min) < fract2) || ((max- u) < fract2))
81 delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
87 delta2 = (x >= mid) ? min - x : max - x;
88 delta = golden * delta2;
90 // update current position:
91 u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
95 // good new point is an improvement!
101 // update control points:
111 // Oh dear, point u is worse than what we have already,
112 // even so it *must* be better than one of our endpoints:
117 if((fu <= fw) || (w == x))
119 // however it is at least second best:
125 else if((fu <= fv) || (v == x) || (v == w))
137 return std::make_pair(x, fx);
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>())))
144 boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
145 return brent_find_minima(f, min, max, digits, m);