165 lines · plain
1// (C) Copyright John Maddock 2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6 7#ifndef BOOST_MATH_TOOLS_MINIMA_HPP8#define BOOST_MATH_TOOLS_MINIMA_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/cstdint.hpp>16#include <boost/math/tools/tuple.hpp>17#include <boost/math/tools/numeric_limits.hpp>18#include <boost/math/tools/precision.hpp>19#include <boost/math/tools/utility.hpp>20#include <boost/math/policies/policy.hpp>21 22namespace boost{ namespace math{ namespace tools{23 24template <class F, class T>25BOOST_MATH_GPU_ENABLED boost::math::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::math::uintmax_t& max_iter)26 noexcept(BOOST_MATH_IS_FLOAT(T) 27 #ifndef BOOST_MATH_HAS_GPU_SUPPORT28 && noexcept(std::declval<F>()(std::declval<T>()))29 #endif30 )31{32 BOOST_MATH_STD_USING33 bits = (boost::math::min)(policies::digits<T, policies::policy<> >() / 2, bits);34 T tolerance = static_cast<T>(ldexp(1.0, 1-bits));35 T x; // minima so far36 T w; // second best point37 T v; // previous value of w38 T u; // most recent evaluation point39 T delta; // The distance moved in the last step40 T delta2; // The distance moved in the step before last41 T fu, fv, fw, fx; // function evaluations at u, v, w, x42 T mid; // midpoint of min and max43 T fract1, fract2; // minimal relative movement in x44 45 static const T golden = 0.3819660f; // golden ratio, don't need too much precision here!46 47 x = w = v = max;48 fw = fv = fx = f(x);49 delta2 = delta = 0;50 51 boost::math::uintmax_t count = max_iter;52 53 do{54 // get midpoint55 mid = (min + max) / 2;56 // work out if we're done already:57 fract1 = tolerance * fabs(x) + tolerance / 4;58 fract2 = 2 * fract1;59 if(fabs(x - mid) <= (fract2 - (max - min) / 2))60 break;61 62 if(fabs(delta2) > fract1)63 {64 // try and construct a parabolic fit:65 T r = (x - w) * (fx - fv);66 T q = (x - v) * (fx - fw);67 T p = (x - v) * q - (x - w) * r;68 q = 2 * (q - r);69 if(q > 0)70 p = -p;71 q = fabs(q);72 T td = delta2;73 delta2 = delta;74 // determine whether a parabolic step is acceptable or not:75 if((fabs(p) >= fabs(q * td / 2)) || (p <= q * (min - x)) || (p >= q * (max - x)))76 {77 // nope, try golden section instead78 delta2 = (x >= mid) ? min - x : max - x;79 delta = golden * delta2;80 }81 else82 {83 // whew, parabolic fit:84 delta = p / q;85 u = x + delta;86 if(((u - min) < fract2) || ((max- u) < fract2))87 delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);88 }89 }90 else91 {92 // golden section:93 delta2 = (x >= mid) ? min - x : max - x;94 delta = golden * delta2;95 }96 // update current position:97 u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));98 fu = f(u);99 if(fu <= fx)100 {101 // good new point is an improvement!102 // update brackets:103 if(u >= x)104 min = x;105 else106 max = x;107 // update control points:108 v = w;109 w = x;110 x = u;111 fv = fw;112 fw = fx;113 fx = fu;114 }115 else116 {117 // Oh dear, point u is worse than what we have already,118 // even so it *must* be better than one of our endpoints:119 if(u < x)120 min = u;121 else122 max = u;123 if((fu <= fw) || (w == x))124 {125 // however it is at least second best:126 v = w;127 w = u;128 fv = fw;129 fw = fu;130 }131 else if((fu <= fv) || (v == x) || (v == w))132 {133 // third best:134 v = u;135 fv = fu;136 }137 }138 139 }while(--count);140 141 max_iter -= count;142 143 return boost::math::make_pair(x, fx);144}145 146template <class F, class T>147BOOST_MATH_GPU_ENABLED inline boost::math::pair<T, T> brent_find_minima(F f, T min, T max, int digits)148 noexcept(BOOST_MATH_IS_FLOAT(T)149 #ifndef BOOST_MATH_HAS_GPU_SUPPORT150 && noexcept(std::declval<F>()(std::declval<T>()))151 #endif152 )153{154 boost::math::uintmax_t m = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();155 return brent_find_minima(f, min, max, digits, m);156}157 158}}} // namespaces159 160#endif161 162 163 164 165