]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/example/brent_minimise_example.cpp
import new upstream nautilus stable release 14.2.8
[ceph.git] / ceph / src / boost / libs / math / example / brent_minimise_example.cpp
CommitLineData
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>
45using 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
62double f(double x)
63{
64 return (x + 3) * (x - 1) * (x - 1);
65}
66
67//[brent_minimise_double_functor
68struct 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
78struct 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.
89template<typename FPT>
90inline bool
91is_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
103template <class T>
104bool 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
124template <class T>
125void 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
187int 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
399For type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>,
400epsilon = 5.3455294202e-51,
401the maximum theoretical precision from Brent minimization is 7.311312755e-26
402Displaying to std::numeric_limits<T>::digits10 11 significant decimal digits.
403x at minimum = 1, f(1) = 5.6273022713e-58,
404met 84 bits precision, after 14 iterations.
405x == 1 (compared to uncertainty 7.311312755e-26) is true
406f(x) == (0 compared to uncertainty 7.311312755e-26) is true
407-4 1.3333333333333333333333333333333333333333333333333
408x at minimum = 0.99999999999999999999999999998813903221565569205253,
409f(0.99999999999999999999999999998813903221565569205253) =
410 5.6273022712501408640665300316078046703496236636624e-58
41114 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
535Typical output MSVC 15.7.3
536
537brent_minimise_example.cpp
538Generating code
5397 of 2746 functions ( 0.3%) were compiled, the rest were copied from previous compilation.
5400 functions were new in current compilation
5411 functions had inline decision re-evaluated but remain unchanged
542Finished generating code
543brent_minimise_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\brent_minimise_example.exe
544Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\brent_minimise_example.exe"
545Brent's minimisation examples.
546
547
548
549Type double - unlimited iterations (unwise?)
550x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-18
551x at minimum = 1.12344622367552e-09
552Uncertainty sqrt(epsilon) = 1.49011611938477e-08
553x = 1.00000, f(x) = 5.04853e-18
554x == 1 (compared to uncertainty 1.49012e-08) is true
555f(x) == 0 (compared to uncertainty 1.49012e-08) is true
556
557Type double with limited iterations.
558Precision bits = 53
559x at minimum = 1.00000, f(1.00000) = 5.04853e-18 after 10 iterations.
560Showing 53 bits precision with 9 decimal digits from tolerance 1.49011612e-08
561x at minimum = 1.00000000, f(1.00000000) = 5.04852568e-18 after 10 iterations.
562
563Type double with limited iterations and half double bits.
564Showing 26 bits precision with 7 decimal digits from tolerance 0.000172633
565x at minimum = 1.000000, f(1.000000) = 5.048526e-18
56610 iterations.
567
568Type double with limited iterations and quarter double bits.
569Showing 13 bits precision with 5 decimal digits from tolerance 0.0156250
570x at minimum = 0.99998, f(0.99998) = 2.0070e-09, after 7 iterations.
571
572Type long double with limited iterations and all long double bits.
573x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-18, after 10 iterations.
574
575
576For type: float,
577epsilon = 1.1921e-07,
578the maximum theoretical precision from Brent's minimization is 0.00034527
579Displaying to std::numeric_limits<T>::digits10 5, significant decimal digits.
580x at minimum = 1.0002, f(1.0002) = 1.9017e-07,
581met 12 bits precision, after 7 iterations.
582x == 1 (compared to uncertainty 0.00034527) is true
583f(x) == (0 compared to uncertainty 0.00034527) is true
584
585
586For type: double,
587epsilon = 2.220446e-16,
588the maximum theoretical precision from Brent's minimization is 1.490116e-08
589Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits.
590x at minimum = 1.000000, f(1.000000) = 5.048526e-18,
591met 26 bits precision, after 10 iterations.
592x == 1 (compared to uncertainty 1.490116e-08) is true
593f(x) == (0 compared to uncertainty 1.490116e-08) is true
594
595
596For type: long double,
597epsilon = 2.220446e-16,
598the maximum theoretical precision from Brent's minimization is 1.490116e-08
599Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits.
600x at minimum = 1.000000, f(1.000000) = 5.048526e-18,
601met 26 bits precision, after 10 iterations.
602x == 1 (compared to uncertainty 1.490116e-08) is true
603f(x) == (0 compared to uncertainty 1.490116e-08) is true
604Bracketing -4.0000000000000000000000000000000000000000000000000 to 1.3333333333333333333333333333333333333333333333333
605x at minimum = 0.99999999999999999999999999998813903221565569205253,
606f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58, after 14 iterations.
607
608
609For type: class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>,
610epsilon = 5.3455294202e-51,
611the maximum theoretical precision from Brent's minimization is 7.3113127550e-26
612Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits.
613x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58,
614met 84 bits precision, after 14 iterations.
615x == 1 (compared to uncertainty 7.3113127550e-26) is true
616f(x) == (0 compared to uncertainty 7.3113127550e-26) is true
617-4.0000000000000000000000000000000000000000000000000 1.3333333333333333333333333333333333333333333333333
618x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58
61914 iterations.
620
7c673cae 621
92f5a8d4
TL
622For type: class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,1>,
623epsilon = 5.3455294202e-51,
624the maximum theoretical precision from Brent's minimization is 7.3113127550e-26
625Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits.
626x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58,
627met 84 bits precision, after 14 iterations.
628x == 1 (compared to uncertainty 7.3113127550e-26) is true
629f(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
636Brent's minimisation examples.
637
638Type double - unlimited iterations (unwise?)
7c673cae
FG
639x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018
640x at minimum = 1.12344622367552e-009
641Uncertainty sqrt(epsilon) = 1.49011611938477e-008
92f5a8d4
TL
642x = 1.00000, f(x) = 5.04853e-018
643x == 1 (compared to uncertainty 1.49012e-008) is true
644f(x) == 0 (compared to uncertainty 1.49012e-008) is true
645
646Type double with limited iterations.
7c673cae 647Precision bits = 53
92f5a8d4
TL
648x at minimum = 1.00000, f(1.00000) = 5.04853e-018 after 10 iterations.
649Showing 53 bits precision with 9 decimal digits from tolerance 1.49011612e-008
650x at minimum = 1.00000000, f(1.00000000) = 5.04852568e-018 after 10 iterations.
651
652Type double with limited iterations and half double bits.
653Showing 26 bits precision with 7 decimal digits from tolerance 0.000172633
654x at minimum = 1.000000, f(1.000000) = 5.048526e-018
7c673cae 65510 iterations.
7c673cae 656
92f5a8d4
TL
657Type double with limited iterations and quarter double bits.
658Showing 13 bits precision with 5 decimal digits from tolerance 0.0156250
659x at minimum = 0.99998, f(0.99998) = 2.0070e-009, after 7 iterations.
660
661Type long double with limited iterations and all long double bits.
662x at minimum = 1.00000000000137302, f(1.00000000000137302) = 7.54079013697311930e-024, after 10 iterations.
7c673cae 663
92f5a8d4
TL
664
665For type: f,
666epsilon = 1.1921e-007,
667the maximum theoretical precision from Brent's minimization is 0.00034527
668Displaying to std::numeric_limits<T>::digits10 5, significant decimal digits.
669x at minimum = 1.0002, f(1.0002) = 1.9017e-007,
670met 12 bits precision, after 7 iterations.
7c673cae
FG
671x == 1 (compared to uncertainty 0.00034527) is true
672f(x) == (0 compared to uncertainty 0.00034527) is true
673
674
92f5a8d4
TL
675For type: d,
676epsilon = 2.220446e-016,
677the maximum theoretical precision from Brent's minimization is 1.490116e-008
678Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits.
679x at minimum = 1.000000, f(1.000000) = 5.048526e-018,
680met 26 bits precision, after 10 iterations.
7c673cae
FG
681x == 1 (compared to uncertainty 1.490116e-008) is true
682f(x) == (0 compared to uncertainty 1.490116e-008) is true
683
684
92f5a8d4
TL
685For type: e,
686epsilon = 1.084202e-019,
687the maximum theoretical precision from Brent's minimization is 3.292723e-010
688Displaying to std::numeric_limits<T>::digits10 7, significant decimal digits.
689x at minimum = 1.000000, f(1.000000) = 7.540790e-024,
690met 32 bits precision, after 10 iterations.
7c673cae
FG
691x == 1 (compared to uncertainty 3.292723e-010) is true
692f(x) == (0 compared to uncertainty 3.292723e-010) is true
693
694
92f5a8d4
TL
695For type: N5boost14multiprecision6numberINS0_8backends16float128_backendELNS0_26expression_template_optionE0EEE,
696epsilon = 1.92592994e-34,
697the maximum theoretical precision from Brent's minimization is 1.38777878e-17
698Displaying to std::numeric_limits<T>::digits10 9, significant decimal digits.
699x at minimum = 1.00000000, f(1.00000000) = 1.48695468e-43,
700met 56 bits precision, after 12 iterations.
7c673cae
FG
701x == 1 (compared to uncertainty 1.38777878e-17) is true
702f(x) == (0 compared to uncertainty 1.38777878e-17) is true
92f5a8d4
TL
703Bracketing -4.0000000000000000000000000000000000000000000000000 to 1.3333333333333333333333333333333333333333333333333
704x at minimum = 0.99999999999999999999999999998813903221565569205253,
705f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58, after 14 iterations.
706
707
708For type: N5boost14multiprecision6numberINS0_8backends13cpp_bin_floatILj50ELNS2_15digit_base_typeE10EviLi0ELi0EEELNS0_26expression_template_optionE1EEE,
709epsilon = 5.3455294202e-51,
710the maximum theoretical precision from Brent's minimization is 7.3113127550e-26
711Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits.
712x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58,
713met 84 bits precision, after 14 iterations.
714x == 1 (compared to uncertainty 7.3113127550e-26) is true
715f(x) == (0 compared to uncertainty 7.3113127550e-26) is true
716-4.0000000000000000000000000000000000000000000000000 1.3333333333333333333333333333333333333333333333333
7c673cae
FG
717x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58
71814 iterations.
719
720
92f5a8d4
TL
721For type: N5boost14multiprecision6numberINS0_8backends13cpp_bin_floatILj50ELNS2_15digit_base_typeE10EviLi0ELi0EEELNS0_26expression_template_optionE1EEE,
722epsilon = 5.3455294202e-51,
723the maximum theoretical precision from Brent's minimization is 7.3113127550e-26
724Displaying to std::numeric_limits<T>::digits10 11, significant decimal digits.
725x at minimum = 1.0000000000, f(1.0000000000) = 5.6273022713e-58,
726met 84 bits precision, after 14 iterations.
727x == 1 (compared to uncertainty 7.3113127550e-26) is true
728f(x) == (0 compared to uncertainty 7.3113127550e-26) is true
7c673cae
FG
729
730*/