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
20namespace boost{ namespace math{ namespace tools{
21
22template <class F, class T>
23std::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 acceptable 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
140template <class F, class T>
141inline 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
155

source code of include/boost/math/tools/minima.hpp