]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | //! \file |
2 | //! \brief Brent_minimise_example.cpp | |
3 | ||
92f5a8d4 | 4 | // Copyright Paul A. Bristow 2015, 2018. |
7c673cae FG |
5 | |
6 | // Use, modification and distribution are subject to the | |
7 | // Boost Software License, Version 1.0. | |
8 | // (See accompanying file LICENSE_1_0.txt | |
9 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | // Note that this file contains Quickbook mark-up as well as code | |
12 | // and comments, don't change any of the special comment mark-ups! | |
13 | ||
14 | // For some diagnostic information: | |
15 | //#define BOOST_MATH_INSTRUMENT | |
16 | // If quadmath float128 is available: | |
17 | //#define BOOST_HAVE_QUADMATH | |
18 | ||
19 | // Example of finding minimum of a function with Brent's method. | |
20 | //[brent_minimise_include_1 | |
21 | #include <boost/math/tools/minima.hpp> | |
22 | //] [/brent_minimise_include_1] | |
23 | ||
24 | #include <boost/math/special_functions/next.hpp> | |
25 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
26 | #include <boost/math/special_functions/pow.hpp> | |
27 | #include <boost/math/constants/constants.hpp> | |
92f5a8d4 | 28 | #include <boost/test/tools/floating_point_comparison.hpp> // For is_close_at)tolerance and is_small |
7c673cae FG |
29 | |
30 | //[brent_minimise_mp_include_0 | |
31 | #include <boost/multiprecision/cpp_dec_float.hpp> // For decimal boost::multiprecision::cpp_dec_float_50. | |
32 | #include <boost/multiprecision/cpp_bin_float.hpp> // For binary boost::multiprecision::cpp_bin_float_50; | |
33 | //] [/brent_minimise_mp_include_0] | |
34 | ||
92f5a8d4 | 35 | //#ifndef _MSC_VER // float128 is not yet supported by Microsoft compiler at 2018. |
7c673cae FG |
36 | #ifdef BOOST_HAVE_QUADMATH // Define only if GCC or Intel, and have quadmath.lib or .dll library available. |
37 | # include <boost/multiprecision/float128.hpp> | |
38 | #endif | |
39 | ||
40 | #include <iostream> | |
41 | // using std::cout; using std::endl; | |
42 | #include <iomanip> | |
43 | // using std::setw; using std::setprecision; | |
44 | #include <limits> | |
45 | using std::numeric_limits; | |
46 | #include <tuple> | |
47 | #include <utility> // pair, make_pair | |
48 | #include <type_traits> | |
49 | #include <typeinfo> | |
50 | ||
51 | //typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>, | |
52 | // boost::multiprecision::et_off> | |
53 | // cpp_dec_float_50_et_off; | |
54 | // | |
55 | // typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>, | |
56 | // boost::multiprecision::et_off> | |
57 | // cpp_bin_float_50_et_off; | |
58 | ||
59 | // http://en.wikipedia.org/wiki/Brent%27s_method Brent's method | |
60 | ||
92f5a8d4 | 61 | // An example of a function for which we want to find a minimum. |
7c673cae FG |
62 | double f(double x) |
63 | { | |
64 | return (x + 3) * (x - 1) * (x - 1); | |
65 | } | |
66 | ||
67 | //[brent_minimise_double_functor | |
68 | struct funcdouble | |
69 | { | |
70 | double operator()(double const& x) | |
92f5a8d4 | 71 | { |
7c673cae FG |
72 | return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2 |
73 | } | |
74 | }; | |
75 | //] [/brent_minimise_double_functor] | |
76 | ||
77 | //[brent_minimise_T_functor | |
78 | struct func | |
79 | { | |
80 | template <class T> | |
81 | T operator()(T const& x) | |
92f5a8d4 TL |
82 | { |
83 | return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2 | |
7c673cae FG |
84 | } |
85 | }; | |
86 | //] [/brent_minimise_T_functor] | |
87 | ||
92f5a8d4 TL |
88 | //! Test if two values are close within a given tolerance. |
89 | template<typename FPT> | |
90 | inline bool | |
91 | is_close_to(FPT left, FPT right, FPT tolerance) | |
92 | { | |
93 | return boost::math::fpc::close_at_tolerance<FPT>(tolerance) (left, right); | |
94 | } | |
95 | ||
7c673cae | 96 | //[brent_minimise_close |
92f5a8d4 TL |
97 | |
98 | //! Compare if value got is close to expected, | |
99 | //! checking first if expected is very small | |
100 | //! (to avoid divide by tiny or zero during comparison) | |
101 | //! before comparing expect with value got. | |
102 | ||
103 | template <class T> | |
104 | bool is_close(T expect, T got, T tolerance) | |
7c673cae | 105 | { |
92f5a8d4 | 106 | using boost::math::fpc::close_at_tolerance; |
7c673cae | 107 | using boost::math::fpc::is_small; |
92f5a8d4 | 108 | using boost::math::fpc::FPC_STRONG; |
7c673cae FG |
109 | |
110 | if (is_small<T>(expect, tolerance)) | |
111 | { | |
112 | return is_small<T>(got, tolerance); | |
113 | } | |
92f5a8d4 TL |
114 | |
115 | return close_at_tolerance<T>(tolerance, FPC_STRONG) (expect, got); | |
116 | } // bool is_close(T expect, T got, T tolerance) | |
7c673cae FG |
117 | |
118 | //] [/brent_minimise_close] | |
119 | ||
120 | //[brent_minimise_T_show | |
121 | ||
92f5a8d4 TL |
122 | //! Example template function to find and show minima. |
123 | //! \tparam T floating-point or fixed_point type. | |
7c673cae FG |
124 | template <class T> |
125 | void show_minima() | |
126 | { | |
127 | using boost::math::tools::brent_find_minima; | |
92f5a8d4 | 128 | using std::sqrt; |
7c673cae | 129 | try |
92f5a8d4 | 130 | { // Always use try'n'catch blocks with Boost.Math to ensure you get any error messages. |
7c673cae FG |
131 | |
132 | int bits = std::numeric_limits<T>::digits/2; // Maximum is digits/2; | |
92f5a8d4 TL |
133 | std::streamsize prec = static_cast<int>(2 + sqrt((double)bits)); // Number of significant decimal digits. |
134 | std::streamsize precision = std::cout.precision(prec); // Save and set. | |
7c673cae | 135 | |
92f5a8d4 | 136 | std::cout << "\n\nFor type: " << typeid(T).name() |
7c673cae FG |
137 | << ",\n epsilon = " << std::numeric_limits<T>::epsilon() |
138 | // << ", precision of " << bits << " bits" | |
92f5a8d4 TL |
139 | << ",\n the maximum theoretical precision from Brent's minimization is " |
140 | << sqrt(std::numeric_limits<T>::epsilon()) | |
141 | << "\n Displaying to std::numeric_limits<T>::digits10 " << prec << ", significant decimal digits." | |
7c673cae FG |
142 | << std::endl; |
143 | ||
144 | const boost::uintmax_t maxit = 20; | |
145 | boost::uintmax_t it = maxit; | |
146 | // Construct using string, not double, avoids loss of precision. | |
147 | //T bracket_min = static_cast<T>("-4"); | |
148 | //T bracket_max = static_cast<T>("1.3333333333333333333333333333333333333333333333333"); | |
149 | ||
92f5a8d4 | 150 | // Construction from double may cause loss of precision for multiprecision types like cpp_bin_float, |
7c673cae FG |
151 | // but brackets values are good enough for using Brent minimization. |
152 | T bracket_min = static_cast<T>(-4); | |
153 | T bracket_max = static_cast<T>(1.3333333333333333333333333333333333333333333333333); | |
154 | ||
155 | std::pair<T, T> r = brent_find_minima<func, T>(func(), bracket_min, bracket_max, bits, it); | |
156 | ||
157 | std::cout << " x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second; | |
158 | if (it < maxit) | |
159 | { | |
160 | std::cout << ",\n met " << bits << " bits precision" << ", after " << it << " iterations." << std::endl; | |
161 | } | |
162 | else | |
163 | { | |
164 | std::cout << ",\n did NOT meet " << bits << " bits precision" << " after " << it << " iterations!" << std::endl; | |
165 | } | |
166 | // Check that result is that expected (compared to theoretical uncertainty). | |
167 | T uncertainty = sqrt(std::numeric_limits<T>::epsilon()); | |
92f5a8d4 TL |
168 | std::cout << std::boolalpha << "x == 1 (compared to uncertainty " << uncertainty << ") is " |
169 | << is_close(static_cast<T>(1), r.first, uncertainty) << std::endl; | |
170 | std::cout << std::boolalpha << "f(x) == (0 compared to uncertainty " << uncertainty << ") is " | |
171 | << is_close(static_cast<T>(0), r.second, uncertainty) << std::endl; | |
7c673cae FG |
172 | // Problems with this using multiprecision with expression template on? |
173 | std::cout.precision(precision); // Restore. | |
174 | } | |
175 | catch (const std::exception& e) | |
176 | { // Always useful to include try & catch blocks because default policies | |
177 | // are to throw exceptions on arguments that cause errors like underflow, overflow. | |
178 | // Lacking try & catch blocks, the program will abort without a message below, | |
179 | // which may give some helpful clues as to the cause of the exception. | |
180 | std::cout << | |
181 | "\n""Message from thrown exception was:\n " << e.what() << std::endl; | |
182 | } | |
183 | } // void show_minima() | |
184 | ||
185 | //] [/brent_minimise_T_show] | |
186 | ||
187 | int main() | |
188 | { | |
92f5a8d4 TL |
189 | using boost::math::tools::brent_find_minima; |
190 | using std::sqrt; | |
191 | std::cout << "Brent's minimisation examples." << std::endl; | |
7c673cae | 192 | std::cout << std::boolalpha << std::endl; |
92f5a8d4 | 193 | std::cout << std::showpoint << std::endl; // Show trailing zeros. |
7c673cae FG |
194 | |
195 | // Tip - using | |
196 | // std::cout.precision(std::numeric_limits<T>::digits10); | |
92f5a8d4 TL |
197 | // during debugging is wise because it warns |
198 | // if construction of multiprecision involves conversion from double | |
199 | // by finding random or zero digits after 17th decimal digit. | |
7c673cae | 200 | |
92f5a8d4 TL |
201 | // Specific type double - unlimited iterations (unwise?). |
202 | { | |
203 | std::cout << "\nType double - unlimited iterations (unwise?)" << std::endl; | |
7c673cae | 204 | //[brent_minimise_double_1 |
92f5a8d4 TL |
205 | const int double_bits = std::numeric_limits<double>::digits; |
206 | std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, double_bits); | |
207 | ||
208 | std::streamsize precision_1 = std::cout.precision(std::numeric_limits<double>::digits10); | |
209 | // Show all double precision decimal digits and trailing zeros. | |
210 | std::cout << "x at minimum = " << r.first | |
211 | << ", f(" << r.first << ") = " << r.second << std::endl; | |
212 | //] [/brent_minimise_double_1] | |
213 | std::cout << "x at minimum = " << (r.first - 1.) / r.first << std::endl; | |
214 | // x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018 | |
215 | double uncertainty = sqrt(std::numeric_limits<double>::epsilon()); | |
216 | std::cout << "Uncertainty sqrt(epsilon) = " << uncertainty << std::endl; | |
217 | // sqrt(epsilon) = 1.49011611938477e-008 | |
218 | // (epsilon is always > 0, so no need to take abs value). | |
219 | ||
220 | std::cout.precision(precision_1); // Restore. | |
221 | //[brent_minimise_double_1a | |
222 | ||
223 | using boost::math::fpc::close_at_tolerance; | |
7c673cae FG |
224 | using boost::math::fpc::is_small; |
225 | ||
92f5a8d4 TL |
226 | std::cout << "x = " << r.first << ", f(x) = " << r.second << std::endl; |
227 | std::cout << std::boolalpha << "x == 1 (compared to uncertainty " | |
228 | << uncertainty << ") is " << is_close(1., r.first, uncertainty) << std::endl; // true | |
229 | std::cout << std::boolalpha << "f(x) == 0 (compared to uncertainty " | |
230 | << uncertainty << ") is " << is_close(0., r.second, uncertainty) << std::endl; // true | |
231 | //] [/brent_minimise_double_1a] | |
7c673cae | 232 | |
92f5a8d4 TL |
233 | } |
234 | std::cout << "\nType double with limited iterations." << std::endl; | |
235 | { | |
236 | const int bits = std::numeric_limits<double>::digits; | |
237 | // Specific type double - limit maxit to 20 iterations. | |
238 | std::cout << "Precision bits = " << bits << std::endl; | |
239 | //[brent_minimise_double_2 | |
240 | const boost::uintmax_t maxit = 20; | |
241 | boost::uintmax_t it = maxit; | |
242 | std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits, it); | |
243 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second | |
244 | << " after " << it << " iterations. " << std::endl; | |
245 | //] [/brent_minimise_double_2] | |
246 | // x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018 | |
7c673cae | 247 | //[brent_minimise_double_3 |
92f5a8d4 TL |
248 | std::streamsize prec = static_cast<int>(2 + sqrt((double)bits)); // Number of significant decimal digits. |
249 | std::streamsize precision_3 = std::cout.precision(prec); // Save and set new precision. | |
250 | std::cout << "Showing " << bits << " bits " | |
251 | "precision with " << prec | |
7c673cae FG |
252 | << " decimal digits from tolerance " << sqrt(std::numeric_limits<double>::epsilon()) |
253 | << std::endl; | |
7c673cae | 254 | |
92f5a8d4 TL |
255 | std::cout << "x at minimum = " << r.first |
256 | << ", f(" << r.first << ") = " << r.second | |
7c673cae | 257 | << " after " << it << " iterations. " << std::endl; |
92f5a8d4 | 258 | std::cout.precision(precision_3); // Restore. |
7c673cae FG |
259 | //] [/brent_minimise_double_3] |
260 | // Showing 53 bits precision with 9 decimal digits from tolerance 1.49011611938477e-008 | |
261 | // x at minimum = 1, f(1) = 5.04852568e-018 | |
92f5a8d4 | 262 | } |
7c673cae | 263 | |
92f5a8d4 | 264 | std::cout << "\nType double with limited iterations and half double bits." << std::endl; |
7c673cae | 265 | { |
92f5a8d4 | 266 | |
7c673cae | 267 | //[brent_minimise_double_4 |
92f5a8d4 | 268 | const int bits_div_2 = std::numeric_limits<double>::digits / 2; // Half digits precision (effective maximum). |
7c673cae | 269 | double epsilon_2 = boost::math::pow<-(std::numeric_limits<double>::digits/2 - 1), double>(2); |
92f5a8d4 | 270 | std::streamsize prec = static_cast<int>(2 + sqrt((double)bits_div_2)); // Number of significant decimal digits. |
7c673cae | 271 | |
92f5a8d4 | 272 | std::cout << "Showing " << bits_div_2 << " bits precision with " << prec |
7c673cae FG |
273 | << " decimal digits from tolerance " << sqrt(epsilon_2) |
274 | << std::endl; | |
92f5a8d4 TL |
275 | std::streamsize precision_4 = std::cout.precision(prec); // Save. |
276 | const boost::uintmax_t maxit = 20; | |
277 | boost::uintmax_t it_4 = maxit; | |
278 | std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits_div_2, it_4); | |
7c673cae | 279 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; |
92f5a8d4 TL |
280 | std::cout << it_4 << " iterations. " << std::endl; |
281 | std::cout.precision(precision_4); // Restore. | |
282 | ||
7c673cae FG |
283 | //] [/brent_minimise_double_4] |
284 | } | |
285 | // x at minimum = 1, f(1) = 5.04852568e-018 | |
286 | ||
287 | { | |
92f5a8d4 | 288 | std::cout << "\nType double with limited iterations and quarter double bits." << std::endl; |
7c673cae | 289 | //[brent_minimise_double_5 |
92f5a8d4 | 290 | const int bits_div_4 = std::numeric_limits<double>::digits / 4; // Quarter precision. |
7c673cae | 291 | double epsilon_4 = boost::math::pow<-(std::numeric_limits<double>::digits / 4 - 1), double>(2); |
92f5a8d4 TL |
292 | std::streamsize prec = static_cast<int>(2 + sqrt((double)bits_div_4)); // Number of significant decimal digits. |
293 | std::cout << "Showing " << bits_div_4 << " bits precision with " << prec | |
7c673cae FG |
294 | << " decimal digits from tolerance " << sqrt(epsilon_4) |
295 | << std::endl; | |
92f5a8d4 TL |
296 | std::streamsize precision_5 = std::cout.precision(prec); // Save & set. |
297 | const boost::uintmax_t maxit = 20; | |
7c673cae | 298 | |
92f5a8d4 TL |
299 | boost::uintmax_t it_5 = maxit; |
300 | std::pair<double, double> r = brent_find_minima(funcdouble(), -4., 4. / 3, bits_div_4, it_5); | |
7c673cae | 301 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second |
92f5a8d4 TL |
302 | << ", after " << it_5 << " iterations. " << std::endl; |
303 | std::cout.precision(precision_5); // Restore. | |
304 | ||
7c673cae FG |
305 | //] [/brent_minimise_double_5] |
306 | } | |
307 | ||
308 | // Showing 13 bits precision with 9 decimal digits from tolerance 0.015625 | |
309 | // x at minimum = 0.9999776, f(0.9999776) = 2.0069572e-009 | |
310 | // 7 iterations. | |
311 | ||
312 | { | |
92f5a8d4 | 313 | std::cout << "\nType long double with limited iterations and all long double bits." << std::endl; |
7c673cae | 314 | //[brent_minimise_template_1 |
92f5a8d4 | 315 | std::streamsize precision_t1 = std::cout.precision(std::numeric_limits<long double>::digits10); // Save & set. |
7c673cae FG |
316 | long double bracket_min = -4.; |
317 | long double bracket_max = 4. / 3; | |
92f5a8d4 | 318 | const int bits = std::numeric_limits<long double>::digits; |
7c673cae FG |
319 | const boost::uintmax_t maxit = 20; |
320 | boost::uintmax_t it = maxit; | |
321 | ||
322 | std::pair<long double, long double> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
323 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second | |
324 | << ", after " << it << " iterations. " << std::endl; | |
92f5a8d4 | 325 | std::cout.precision(precision_t1); // Restore. |
7c673cae FG |
326 | //] [/brent_minimise_template_1] |
327 | } | |
328 | ||
329 | // Show use of built-in type Template versions. | |
330 | // (Will not work if construct bracket min and max from string). | |
331 | ||
332 | //[brent_minimise_template_fd | |
333 | show_minima<float>(); | |
334 | show_minima<double>(); | |
335 | show_minima<long double>(); | |
7c673cae | 336 | |
92f5a8d4 | 337 | //] [/brent_minimise_template_fd] |
7c673cae FG |
338 | |
339 | //[brent_minimise_mp_include_1 | |
92f5a8d4 | 340 | #ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available. |
7c673cae FG |
341 | using boost::multiprecision::float128; |
342 | #endif | |
343 | //] [/brent_minimise_mp_include_1] | |
344 | ||
345 | //[brent_minimise_template_quad | |
92f5a8d4 | 346 | #ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available. |
7c673cae FG |
347 | show_minima<float128>(); // Needs quadmath_snprintf, sqrtQ, fabsq that are in in quadmath library. |
348 | #endif | |
349 | //] [/brent_minimise_template_quad | |
350 | ||
351 | // User-defined floating-point template. | |
352 | ||
353 | //[brent_minimise_mp_typedefs | |
92f5a8d4 TL |
354 | using boost::multiprecision::cpp_bin_float_50; // binary multiprecision typedef. |
355 | using boost::multiprecision::cpp_dec_float_50; // decimal multiprecision typedef. | |
7c673cae | 356 | |
92f5a8d4 | 357 | // One might also need typedefs like these to switch expression templates off and on (default is on). |
7c673cae FG |
358 | typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>, |
359 | boost::multiprecision::et_on> | |
360 | cpp_bin_float_50_et_on; // et_on is default so is same as cpp_bin_float_50. | |
361 | ||
362 | typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>, | |
363 | boost::multiprecision::et_off> | |
364 | cpp_bin_float_50_et_off; | |
365 | ||
7c673cae FG |
366 | typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>, |
367 | boost::multiprecision::et_on> // et_on is default so is same as cpp_dec_float_50. | |
368 | cpp_dec_float_50_et_on; | |
369 | ||
370 | typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>, | |
371 | boost::multiprecision::et_off> | |
372 | cpp_dec_float_50_et_off; | |
373 | //] [/brent_minimise_mp_typedefs] | |
374 | ||
375 | { // binary ET on by default. | |
376 | //[brent_minimise_mp_1 | |
377 | std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10); | |
7c673cae | 378 | int bits = std::numeric_limits<cpp_bin_float_50>::digits / 2 - 2; |
7c673cae FG |
379 | cpp_bin_float_50 bracket_min = static_cast<cpp_bin_float_50>("-4"); |
380 | cpp_bin_float_50 bracket_max = static_cast<cpp_bin_float_50>("1.3333333333333333333333333333333333333333333333333"); | |
381 | ||
92f5a8d4 | 382 | std::cout << "Bracketing " << bracket_min << " to " << bracket_max << std::endl; |
7c673cae | 383 | const boost::uintmax_t maxit = 20; |
92f5a8d4 TL |
384 | boost::uintmax_t it = maxit; // Will be updated with actual iteration count. |
385 | std::pair<cpp_bin_float_50, cpp_bin_float_50> r | |
386 | = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
7c673cae | 387 | |
92f5a8d4 | 388 | std::cout << "x at minimum = " << r.first << ",\n f(" << r.first << ") = " << r.second |
7c673cae FG |
389 | // x at minimum = 1, f(1) = 5.04853e-018 |
390 | << ", after " << it << " iterations. " << std::endl; | |
391 | ||
92f5a8d4 TL |
392 | is_close_to(static_cast<cpp_bin_float_50>("1"), r.first, sqrt(std::numeric_limits<cpp_bin_float_50>::epsilon())); |
393 | is_close_to(static_cast<cpp_bin_float_50>("0"), r.second, sqrt(std::numeric_limits<cpp_bin_float_50>::epsilon())); | |
7c673cae FG |
394 | |
395 | //] [/brent_minimise_mp_1] | |
396 | ||
397 | /* | |
398 | //[brent_minimise_mp_output_1 | |
92f5a8d4 TL |
399 | For type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>, |
400 | epsilon = 5.3455294202e-51, | |
401 | the maximum theoretical precision from Brent minimization is 7.311312755e-26 | |
402 | Displaying to std::numeric_limits<T>::digits10 11 significant decimal digits. | |
403 | x at minimum = 1, f(1) = 5.6273022713e-58, | |
404 | met 84 bits precision, after 14 iterations. | |
405 | x == 1 (compared to uncertainty 7.311312755e-26) is true | |
406 | f(x) == (0 compared to uncertainty 7.311312755e-26) is true | |
407 | -4 1.3333333333333333333333333333333333333333333333333 | |
408 | x at minimum = 0.99999999999999999999999999998813903221565569205253, | |
409 | f(0.99999999999999999999999999998813903221565569205253) = | |
410 | 5.6273022712501408640665300316078046703496236636624e-58 | |
411 | 14 iterations | |
7c673cae FG |
412 | //] [/brent_minimise_mp_output_1] |
413 | */ | |
414 | //[brent_minimise_mp_2 | |
415 | show_minima<cpp_bin_float_50_et_on>(); // | |
416 | //] [/brent_minimise_mp_2] | |
417 | ||
418 | /* | |
419 | //[brent_minimise_mp_output_2 | |
420 | For type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50, 10, void, int, 0, 0>, 1>, | |
421 | ||
422 | //] [/brent_minimise_mp_output_1] | |
423 | */ | |
424 | } | |
425 | ||
426 | { // binary ET on explicit | |
427 | std::cout.precision(std::numeric_limits<cpp_bin_float_50_et_on>::digits10); | |
428 | ||
429 | int bits = std::numeric_limits<cpp_bin_float_50_et_on>::digits / 2 - 2; | |
430 | ||
431 | cpp_bin_float_50_et_on bracket_min = static_cast<cpp_bin_float_50_et_on>("-4"); | |
432 | cpp_bin_float_50_et_on bracket_max = static_cast<cpp_bin_float_50_et_on>("1.3333333333333333333333333333333333333333333333333"); | |
433 | ||
434 | std::cout << bracket_min << " " << bracket_max << std::endl; | |
435 | const boost::uintmax_t maxit = 20; | |
436 | boost::uintmax_t it = maxit; | |
437 | std::pair<cpp_bin_float_50_et_on, cpp_bin_float_50_et_on> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
438 | ||
439 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; | |
440 | // x at minimum = 1, f(1) = 5.04853e-018 | |
441 | std::cout << it << " iterations. " << std::endl; | |
442 | ||
443 | show_minima<cpp_bin_float_50_et_on>(); // | |
444 | ||
445 | } | |
446 | return 0; | |
447 | ||
92f5a8d4 TL |
448 | // Some examples of switching expression templates on and off follow. |
449 | ||
7c673cae FG |
450 | { // binary ET off |
451 | std::cout.precision(std::numeric_limits<cpp_bin_float_50_et_off>::digits10); | |
452 | ||
453 | int bits = std::numeric_limits<cpp_bin_float_50_et_off>::digits / 2 - 2; | |
454 | cpp_bin_float_50_et_off bracket_min = static_cast<cpp_bin_float_50_et_off>("-4"); | |
455 | cpp_bin_float_50_et_off bracket_max = static_cast<cpp_bin_float_50_et_off>("1.3333333333333333333333333333333333333333333333333"); | |
456 | ||
457 | std::cout << bracket_min << " " << bracket_max << std::endl; | |
458 | const boost::uintmax_t maxit = 20; | |
459 | boost::uintmax_t it = maxit; | |
460 | std::pair<cpp_bin_float_50_et_off, cpp_bin_float_50_et_off> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
461 | ||
462 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; | |
463 | // x at minimum = 1, f(1) = 5.04853e-018 | |
464 | std::cout << it << " iterations. " << std::endl; | |
465 | ||
466 | show_minima<cpp_bin_float_50_et_off>(); // | |
467 | } | |
468 | ||
469 | { // decimal ET on by default | |
470 | std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10); | |
471 | ||
472 | int bits = std::numeric_limits<cpp_dec_float_50>::digits / 2 - 2; | |
473 | ||
474 | cpp_dec_float_50 bracket_min = static_cast<cpp_dec_float_50>("-4"); | |
475 | cpp_dec_float_50 bracket_max = static_cast<cpp_dec_float_50>("1.3333333333333333333333333333333333333333333333333"); | |
476 | ||
477 | std::cout << bracket_min << " " << bracket_max << std::endl; | |
478 | const boost::uintmax_t maxit = 20; | |
479 | boost::uintmax_t it = maxit; | |
480 | std::pair<cpp_dec_float_50, cpp_dec_float_50> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
481 | ||
482 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; | |
483 | // x at minimum = 1, f(1) = 5.04853e-018 | |
484 | std::cout << it << " iterations. " << std::endl; | |
485 | ||
486 | show_minima<cpp_dec_float_50>(); | |
487 | } | |
488 | ||
489 | { // decimal ET on | |
490 | std::cout.precision(std::numeric_limits<cpp_dec_float_50_et_on>::digits10); | |
491 | ||
492 | int bits = std::numeric_limits<cpp_dec_float_50_et_on>::digits / 2 - 2; | |
493 | ||
494 | cpp_dec_float_50_et_on bracket_min = static_cast<cpp_dec_float_50_et_on>("-4"); | |
495 | cpp_dec_float_50_et_on bracket_max = static_cast<cpp_dec_float_50_et_on>("1.3333333333333333333333333333333333333333333333333"); | |
496 | std::cout << bracket_min << " " << bracket_max << std::endl; | |
497 | const boost::uintmax_t maxit = 20; | |
498 | boost::uintmax_t it = maxit; | |
499 | std::pair<cpp_dec_float_50_et_on, cpp_dec_float_50_et_on> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
500 | ||
501 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; | |
502 | // x at minimum = 1, f(1) = 5.04853e-018 | |
503 | std::cout << it << " iterations. " << std::endl; | |
504 | ||
505 | show_minima<cpp_dec_float_50_et_on>(); | |
506 | ||
507 | } | |
508 | ||
509 | { // decimal ET off | |
510 | std::cout.precision(std::numeric_limits<cpp_dec_float_50_et_off>::digits10); | |
511 | ||
512 | int bits = std::numeric_limits<cpp_dec_float_50_et_off>::digits / 2 - 2; | |
513 | ||
514 | cpp_dec_float_50_et_off bracket_min = static_cast<cpp_dec_float_50_et_off>("-4"); | |
515 | cpp_dec_float_50_et_off bracket_max = static_cast<cpp_dec_float_50_et_off>("1.3333333333333333333333333333333333333333333333333"); | |
516 | ||
517 | std::cout << bracket_min << " " << bracket_max << std::endl; | |
518 | const boost::uintmax_t maxit = 20; | |
519 | boost::uintmax_t it = maxit; | |
520 | std::pair<cpp_dec_float_50_et_off, cpp_dec_float_50_et_off> r = brent_find_minima(func(), bracket_min, bracket_max, bits, it); | |
521 | ||
522 | std::cout << "x at minimum = " << r.first << ", f(" << r.first << ") = " << r.second << std::endl; | |
523 | // x at minimum = 1, f(1) = 5.04853e-018 | |
524 | std::cout << it << " iterations. " << std::endl; | |
525 | ||
526 | show_minima<cpp_dec_float_50_et_off>(); | |
527 | } | |
528 | ||
529 | return 0; | |
530 | } // int main() | |
531 | ||
532 | ||
533 | /* | |
534 | ||
92f5a8d4 TL |
535 | Typical output MSVC 15.7.3 |
536 | ||
537 | brent_minimise_example.cpp | |
538 | Generating code | |
539 | 7 of 2746 functions ( 0.3%) were compiled, the rest were copied from previous compilation. | |
540 | 0 functions were new in current compilation | |
541 | 1 functions had inline decision re-evaluated but remain unchanged | |
542 | Finished generating code | |
543 | brent_minimise_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\brent_minimise_example.exe | |
544 | Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\brent_minimise_example.exe" | |
545 | Brent's minimisation examples. | |
546 | ||
547 | ||
548 | ||
549 | Type double - unlimited iterations (unwise?) | |
550 | x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-18 | |
551 | x at minimum = 1.12344622367552e-09 | |
552 | Uncertainty sqrt(epsilon) = 1.49011611938477e-08 | |
553 | x = 1.00000, f(x) = 5.04853e-18 | |
554 | x == 1 (compared to uncertainty 1.49012e-08) is true | |
555 | f(x) == 0 (compared to uncertainty 1.49012e-08) is true | |
556 | ||
557 | Type double with limited iterations. | |
558 | Precision bits = 53 | |
559 | x at minimum = 1.00000, f(1.00000) = 5.04853e-18 after 10 iterations. | |
560 | Showing 53 bits precision with 9 decimal digits from tolerance 1.49011612e-08 | |
561 | x at minimum = 1.00000000, f(1.00000000) = 5.04852568e-18 after 10 iterations. | |
562 | ||
563 | Type double with limited iterations and half double bits. | |
564 | Showing 26 bits precision with 7 decimal digits from tolerance 0.000172633 | |
565 | x at minimum = 1.000000, f(1.000000) = 5.048526e-18 | |
566 | 10 iterations. | |
567 | ||
568 | Type double with limited iterations and quarter double bits. | |
569 | Showing 13 bits precision with 5 decimal digits from tolerance 0.0156250 | |
570 | x at minimum = 0.99998, f(0.99998) = 2.0070e-09, after 7 iterations. | |
571 | ||
572 | Type long double with limited iterations and all long double bits. | |
573 | x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-18, after 10 iterations. | |
574 | ||
575 | ||
576 | For type: float, | |
577 | epsilon = 1.1921e-07, | |
578 | the maximum theoretical precision from Brent's minimization is 0.00034527 | |
579 | Displaying to std::numeric_limits<T>::digits10 5, significant decimal digits. | |
580 | x at minimum = 1.0002, f(1.0002) = 1.9017e-07, | |
581 | met 12 bits precision, after 7 iterations. | |
582 | x == 1 (compared to uncertainty 0.00034527) is true | |
583 | f(x) == (0 compared to uncertainty 0.00034527) is true | |
584 | ||
585 | ||
586 | For type: double, | |
587 | epsilon = 2.220446e-16, | |
588 | the maximum theoretical precision from Brent's minimization is 1.490116e-08 | |
589 | Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits. | |
590 | x at minimum = 1.000000, f(1.000000) = 5.048526e-18, | |
591 | met 26 bits precision, after 10 iterations. | |
592 | x == 1 (compared to uncertainty 1.490116e-08) is true | |
593 | f(x) == (0 compared to uncertainty 1.490116e-08) is true | |
594 | ||
595 | ||
596 | For type: long double, | |
597 | epsilon = 2.220446e-16, | |
598 | the maximum theoretical precision from Brent's minimization is 1.490116e-08 | |
599 | Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits. | |
600 | x at minimum = 1.000000, f(1.000000) = 5.048526e-18, | |
601 | met 26 bits precision, after 10 iterations. | |
602 | x == 1 (compared to uncertainty 1.490116e-08) is true | |
603 | f(x) == (0 compared to uncertainty 1.490116e-08) is true | |
604 | Bracketing -4.0000000000000000000000000000000000000000000000000 to 1.3333333333333333333333333333333333333333333333333 | |
605 | x at minimum = 0.99999999999999999999999999998813903221565569205253, | |
606 | f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58, after 14 iterations. | |
607 | ||
608 | ||
609 | For type: class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>, | |
610 | epsilon = 5.3455294202e-51, | |
611 | the maximum theoretical precision from Brent's minimization is 7.3113127550e-26 | |
612 | Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits. | |
613 | x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58, | |
614 | met 84 bits precision, after 14 iterations. | |
615 | x == 1 (compared to uncertainty 7.3113127550e-26) is true | |
616 | f(x) == (0 compared to uncertainty 7.3113127550e-26) is true | |
617 | -4.0000000000000000000000000000000000000000000000000 1.3333333333333333333333333333333333333333333333333 | |
618 | x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58 | |
619 | 14 iterations. | |
620 | ||
7c673cae | 621 | |
92f5a8d4 TL |
622 | For type: class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>, |
623 | epsilon = 5.3455294202e-51, | |
624 | the maximum theoretical precision from Brent's minimization is 7.3113127550e-26 | |
625 | Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits. | |
626 | x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58, | |
627 | met 84 bits precision, after 14 iterations. | |
628 | x == 1 (compared to uncertainty 7.3113127550e-26) is true | |
629 | f(x) == (0 compared to uncertainty 7.3113127550e-26) is true | |
7c673cae | 630 | |
7c673cae | 631 | |
92f5a8d4 TL |
632 | ============================================================================================================ |
633 | ||
634 | // GCC 7.2.0 with quadmath | |
635 | ||
636 | Brent's minimisation examples. | |
637 | ||
638 | Type double - unlimited iterations (unwise?) | |
7c673cae FG |
639 | x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018 |
640 | x at minimum = 1.12344622367552e-009 | |
641 | Uncertainty sqrt(epsilon) = 1.49011611938477e-008 | |
92f5a8d4 TL |
642 | x = 1.00000, f(x) = 5.04853e-018 |
643 | x == 1 (compared to uncertainty 1.49012e-008) is true | |
644 | f(x) == 0 (compared to uncertainty 1.49012e-008) is true | |
645 | ||
646 | Type double with limited iterations. | |
7c673cae | 647 | Precision bits = 53 |
92f5a8d4 TL |
648 | x at minimum = 1.00000, f(1.00000) = 5.04853e-018 after 10 iterations. |
649 | Showing 53 bits precision with 9 decimal digits from tolerance 1.49011612e-008 | |
650 | x at minimum = 1.00000000, f(1.00000000) = 5.04852568e-018 after 10 iterations. | |
651 | ||
652 | Type double with limited iterations and half double bits. | |
653 | Showing 26 bits precision with 7 decimal digits from tolerance 0.000172633 | |
654 | x at minimum = 1.000000, f(1.000000) = 5.048526e-018 | |
7c673cae | 655 | 10 iterations. |
7c673cae | 656 | |
92f5a8d4 TL |
657 | Type double with limited iterations and quarter double bits. |
658 | Showing 13 bits precision with 5 decimal digits from tolerance 0.0156250 | |
659 | x at minimum = 0.99998, f(0.99998) = 2.0070e-009, after 7 iterations. | |
660 | ||
661 | Type long double with limited iterations and all long double bits. | |
662 | x at minimum = 1.00000000000137302, f(1.00000000000137302) = 7.54079013697311930e-024, after 10 iterations. | |
7c673cae | 663 | |
92f5a8d4 TL |
664 | |
665 | For type: f, | |
666 | epsilon = 1.1921e-007, | |
667 | the maximum theoretical precision from Brent's minimization is 0.00034527 | |
668 | Displaying to std::numeric_limits<T>::digits10 5, significant decimal digits. | |
669 | x at minimum = 1.0002, f(1.0002) = 1.9017e-007, | |
670 | met 12 bits precision, after 7 iterations. | |
7c673cae FG |
671 | x == 1 (compared to uncertainty 0.00034527) is true |
672 | f(x) == (0 compared to uncertainty 0.00034527) is true | |
673 | ||
674 | ||
92f5a8d4 TL |
675 | For type: d, |
676 | epsilon = 2.220446e-016, | |
677 | the maximum theoretical precision from Brent's minimization is 1.490116e-008 | |
678 | Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits. | |
679 | x at minimum = 1.000000, f(1.000000) = 5.048526e-018, | |
680 | met 26 bits precision, after 10 iterations. | |
7c673cae FG |
681 | x == 1 (compared to uncertainty 1.490116e-008) is true |
682 | f(x) == (0 compared to uncertainty 1.490116e-008) is true | |
683 | ||
684 | ||
92f5a8d4 TL |
685 | For type: e, |
686 | epsilon = 1.084202e-019, | |
687 | the maximum theoretical precision from Brent's minimization is 3.292723e-010 | |
688 | Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits. | |
689 | x at minimum = 1.000000, f(1.000000) = 7.540790e-024, | |
690 | met 32 bits precision, after 10 iterations. | |
7c673cae FG |
691 | x == 1 (compared to uncertainty 3.292723e-010) is true |
692 | f(x) == (0 compared to uncertainty 3.292723e-010) is true | |
693 | ||
694 | ||
92f5a8d4 TL |
695 | For type: N5boost14multiprecision6numberINS0_8backends16float128_backendELNS0_26expression_template_optionE0EEE, |
696 | epsilon = 1.92592994e-34, | |
697 | the maximum theoretical precision from Brent's minimization is 1.38777878e-17 | |
698 | Displaying to std::numeric_limits<T>::digits10 9, significant decimal digits. | |
699 | x at minimum = 1.00000000, f(1.00000000) = 1.48695468e-43, | |
700 | met 56 bits precision, after 12 iterations. | |
7c673cae FG |
701 | x == 1 (compared to uncertainty 1.38777878e-17) is true |
702 | f(x) == (0 compared to uncertainty 1.38777878e-17) is true | |
92f5a8d4 TL |
703 | Bracketing -4.0000000000000000000000000000000000000000000000000 to 1.3333333333333333333333333333333333333333333333333 |
704 | x at minimum = 0.99999999999999999999999999998813903221565569205253, | |
705 | f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58, after 14 iterations. | |
706 | ||
707 | ||
708 | For type: N5boost14multiprecision6numberINS0_8backends13cpp_bin_floatILj50ELNS2_15digit_base_typeE10EviLi0ELi0EEELNS0_26expression_template_optionE1EEE, | |
709 | epsilon = 5.3455294202e-51, | |
710 | the maximum theoretical precision from Brent's minimization is 7.3113127550e-26 | |
711 | Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits. | |
712 | x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58, | |
713 | met 84 bits precision, after 14 iterations. | |
714 | x == 1 (compared to uncertainty 7.3113127550e-26) is true | |
715 | f(x) == (0 compared to uncertainty 7.3113127550e-26) is true | |
716 | -4.0000000000000000000000000000000000000000000000000 1.3333333333333333333333333333333333333333333333333 | |
7c673cae FG |
717 | x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58 |
718 | 14 iterations. | |
719 | ||
720 | ||
92f5a8d4 TL |
721 | For type: N5boost14multiprecision6numberINS0_8backends13cpp_bin_floatILj50ELNS2_15digit_base_typeE10EviLi0ELi0EEELNS0_26expression_template_optionE1EEE, |
722 | epsilon = 5.3455294202e-51, | |
723 | the maximum theoretical precision from Brent's minimization is 7.3113127550e-26 | |
724 | Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits. | |
725 | x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58, | |
726 | met 84 bits precision, after 14 iterations. | |
727 | x == 1 (compared to uncertainty 7.3113127550e-26) is true | |
728 | f(x) == (0 compared to uncertainty 7.3113127550e-26) is true | |
7c673cae FG |
729 | |
730 | */ |