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