2 //! \brief Brent_minimise_example.cpp
4 // Copyright Paul A. Bristow 2015, 2018.
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)
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!
14 // For some diagnostic information:
15 //#define BOOST_MATH_INSTRUMENT
16 // If quadmath float128 is available:
17 //#define BOOST_HAVE_QUADMATH
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]
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>
28 #include <boost/test/tools/floating_point_comparison.hpp> // For is_close_at)tolerance and is_small
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]
35 //#ifndef _MSC_VER // float128 is not yet supported by Microsoft compiler at 2018.
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>
41 // using std::cout; using std::endl;
43 // using std::setw; using std::setprecision;
45 using std::numeric_limits
;
47 #include <utility> // pair, make_pair
48 #include <type_traits>
51 //typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<50>,
52 // boost::multiprecision::et_off>
53 // cpp_dec_float_50_et_off;
55 // typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<50>,
56 // boost::multiprecision::et_off>
57 // cpp_bin_float_50_et_off;
59 // http://en.wikipedia.org/wiki/Brent%27s_method Brent's method
61 // An example of a function for which we want to find a minimum.
64 return (x
+ 3) * (x
- 1) * (x
- 1);
67 //[brent_minimise_double_functor
70 double operator()(double const& x
)
72 return (x
+ 3) * (x
- 1) * (x
- 1); // (x + 3)(x - 1)^2
75 //] [/brent_minimise_double_functor]
77 //[brent_minimise_T_functor
81 T
operator()(T
const& x
)
83 return (x
+ 3) * (x
- 1) * (x
- 1); // (x + 3)(x - 1)^2
86 //] [/brent_minimise_T_functor]
88 //! Test if two values are close within a given tolerance.
89 template<typename FPT
>
91 is_close_to(FPT left
, FPT right
, FPT tolerance
)
93 return boost::math::fpc::close_at_tolerance
<FPT
>(tolerance
) (left
, right
);
96 //[brent_minimise_close
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.
104 bool is_close(T expect
, T got
, T tolerance
)
106 using boost::math::fpc::close_at_tolerance
;
107 using boost::math::fpc::is_small
;
108 using boost::math::fpc::FPC_STRONG
;
110 if (is_small
<T
>(expect
, tolerance
))
112 return is_small
<T
>(got
, tolerance
);
115 return close_at_tolerance
<T
>(tolerance
, FPC_STRONG
) (expect
, got
);
116 } // bool is_close(T expect, T got, T tolerance)
118 //] [/brent_minimise_close]
120 //[brent_minimise_T_show
122 //! Example template function to find and show minima.
123 //! \tparam T floating-point or fixed_point type.
127 using boost::math::tools::brent_find_minima
;
130 { // Always use try'n'catch blocks with Boost.Math to ensure you get any error messages.
132 int bits
= std::numeric_limits
<T
>::digits
/2; // Maximum is digits/2;
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.
136 std::cout
<< "\n\nFor type: " << typeid(T
).name()
137 << ",\n epsilon = " << std::numeric_limits
<T
>::epsilon()
138 // << ", precision of " << bits << " bits"
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."
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");
150 // Construction from double may cause loss of precision for multiprecision types like cpp_bin_float,
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);
155 std::pair
<T
, T
> r
= brent_find_minima
<func
, T
>(func(), bracket_min
, bracket_max
, bits
, it
);
157 std::cout
<< " x at minimum = " << r
.first
<< ", f(" << r
.first
<< ") = " << r
.second
;
160 std::cout
<< ",\n met " << bits
<< " bits precision" << ", after " << it
<< " iterations." << std::endl
;
164 std::cout
<< ",\n did NOT meet " << bits
<< " bits precision" << " after " << it
<< " iterations!" << std::endl
;
166 // Check that result is that expected (compared to theoretical uncertainty).
167 T uncertainty
= sqrt(std::numeric_limits
<T
>::epsilon());
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
;
172 // Problems with this using multiprecision with expression template on?
173 std::cout
.precision(precision
); // Restore.
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.
181 "\n""Message from thrown exception was:\n " << e
.what() << std::endl
;
183 } // void show_minima()
185 //] [/brent_minimise_T_show]
189 using boost::math::tools::brent_find_minima
;
191 std::cout
<< "Brent's minimisation examples." << std::endl
;
192 std::cout
<< std::boolalpha
<< std::endl
;
193 std::cout
<< std::showpoint
<< std::endl
; // Show trailing zeros.
196 // std::cout.precision(std::numeric_limits<T>::digits10);
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.
201 // Specific type double - unlimited iterations (unwise?).
203 std::cout
<< "\nType double - unlimited iterations (unwise?)" << std::endl
;
204 //[brent_minimise_double_1
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
);
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).
220 std::cout
.precision(precision_1
); // Restore.
221 //[brent_minimise_double_1a
223 using boost::math::fpc::close_at_tolerance
;
224 using boost::math::fpc::is_small
;
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]
234 std::cout
<< "\nType double with limited iterations." << std::endl
;
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
247 //[brent_minimise_double_3
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
252 << " decimal digits from tolerance " << sqrt(std::numeric_limits
<double>::epsilon())
255 std::cout
<< "x at minimum = " << r
.first
256 << ", f(" << r
.first
<< ") = " << r
.second
257 << " after " << it
<< " iterations. " << std::endl
;
258 std::cout
.precision(precision_3
); // Restore.
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
264 std::cout
<< "\nType double with limited iterations and half double bits." << std::endl
;
267 //[brent_minimise_double_4
268 const int bits_div_2
= std::numeric_limits
<double>::digits
/ 2; // Half digits precision (effective maximum).
269 double epsilon_2
= boost::math::pow
<-(std::numeric_limits
<double>::digits
/2 - 1), double>(2);
270 std::streamsize prec
= static_cast<int>(2 + sqrt((double)bits_div_2
)); // Number of significant decimal digits.
272 std::cout
<< "Showing " << bits_div_2
<< " bits precision with " << prec
273 << " decimal digits from tolerance " << sqrt(epsilon_2
)
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
);
279 std::cout
<< "x at minimum = " << r
.first
<< ", f(" << r
.first
<< ") = " << r
.second
<< std::endl
;
280 std::cout
<< it_4
<< " iterations. " << std::endl
;
281 std::cout
.precision(precision_4
); // Restore.
283 //] [/brent_minimise_double_4]
285 // x at minimum = 1, f(1) = 5.04852568e-018
288 std::cout
<< "\nType double with limited iterations and quarter double bits." << std::endl
;
289 //[brent_minimise_double_5
290 const int bits_div_4
= std::numeric_limits
<double>::digits
/ 4; // Quarter precision.
291 double epsilon_4
= boost::math::pow
<-(std::numeric_limits
<double>::digits
/ 4 - 1), double>(2);
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
294 << " decimal digits from tolerance " << sqrt(epsilon_4
)
296 std::streamsize precision_5
= std::cout
.precision(prec
); // Save & set.
297 const boost::uintmax_t maxit
= 20;
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
);
301 std::cout
<< "x at minimum = " << r
.first
<< ", f(" << r
.first
<< ") = " << r
.second
302 << ", after " << it_5
<< " iterations. " << std::endl
;
303 std::cout
.precision(precision_5
); // Restore.
305 //] [/brent_minimise_double_5]
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
313 std::cout
<< "\nType long double with limited iterations and all long double bits." << std::endl
;
314 //[brent_minimise_template_1
315 std::streamsize precision_t1
= std::cout
.precision(std::numeric_limits
<long double>::digits10
); // Save & set.
316 long double bracket_min
= -4.;
317 long double bracket_max
= 4. / 3;
318 const int bits
= std::numeric_limits
<long double>::digits
;
319 const boost::uintmax_t maxit
= 20;
320 boost::uintmax_t it
= maxit
;
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
;
325 std::cout
.precision(precision_t1
); // Restore.
326 //] [/brent_minimise_template_1]
329 // Show use of built-in type Template versions.
330 // (Will not work if construct bracket min and max from string).
332 //[brent_minimise_template_fd
333 show_minima
<float>();
334 show_minima
<double>();
335 show_minima
<long double>();
337 //] [/brent_minimise_template_fd]
339 //[brent_minimise_mp_include_1
340 #ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available.
341 using boost::multiprecision::float128
;
343 //] [/brent_minimise_mp_include_1]
345 //[brent_minimise_template_quad
346 #ifdef BOOST_HAVE_QUADMATH // Defined only if GCC or Intel and have quadmath.lib or .dll library available.
347 show_minima
<float128
>(); // Needs quadmath_snprintf, sqrtQ, fabsq that are in in quadmath library.
349 //] [/brent_minimise_template_quad
351 // User-defined floating-point template.
353 //[brent_minimise_mp_typedefs
354 using boost::multiprecision::cpp_bin_float_50
; // binary multiprecision typedef.
355 using boost::multiprecision::cpp_dec_float_50
; // decimal multiprecision typedef.
357 // One might also need typedefs like these to switch expression templates off and on (default is on).
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.
362 typedef boost::multiprecision::number
<boost::multiprecision::cpp_bin_float
<50>,
363 boost::multiprecision::et_off
>
364 cpp_bin_float_50_et_off
;
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
;
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]
375 { // binary ET on by default.
376 //[brent_minimise_mp_1
377 std::cout
.precision(std::numeric_limits
<cpp_bin_float_50
>::digits10
);
378 int bits
= std::numeric_limits
<cpp_bin_float_50
>::digits
/ 2 - 2;
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");
382 std::cout
<< "Bracketing " << bracket_min
<< " to " << bracket_max
<< std::endl
;
383 const boost::uintmax_t maxit
= 20;
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
);
388 std::cout
<< "x at minimum = " << r
.first
<< ",\n f(" << r
.first
<< ") = " << r
.second
389 // x at minimum = 1, f(1) = 5.04853e-018
390 << ", after " << it
<< " iterations. " << std::endl
;
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()));
395 //] [/brent_minimise_mp_1]
398 //[brent_minimise_mp_output_1
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
412 //] [/brent_minimise_mp_output_1]
414 //[brent_minimise_mp_2
415 show_minima
<cpp_bin_float_50_et_on
>(); //
416 //] [/brent_minimise_mp_2]
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>,
422 //] [/brent_minimise_mp_output_1]
426 { // binary ET on explicit
427 std::cout
.precision(std::numeric_limits
<cpp_bin_float_50_et_on
>::digits10
);
429 int bits
= std::numeric_limits
<cpp_bin_float_50_et_on
>::digits
/ 2 - 2;
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");
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
);
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
;
443 show_minima
<cpp_bin_float_50_et_on
>(); //
448 // Some examples of switching expression templates on and off follow.
451 std::cout
.precision(std::numeric_limits
<cpp_bin_float_50_et_off
>::digits10
);
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");
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
);
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
;
466 show_minima
<cpp_bin_float_50_et_off
>(); //
469 { // decimal ET on by default
470 std::cout
.precision(std::numeric_limits
<cpp_dec_float_50
>::digits10
);
472 int bits
= std::numeric_limits
<cpp_dec_float_50
>::digits
/ 2 - 2;
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");
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
);
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
;
486 show_minima
<cpp_dec_float_50
>();
490 std::cout
.precision(std::numeric_limits
<cpp_dec_float_50_et_on
>::digits10
);
492 int bits
= std::numeric_limits
<cpp_dec_float_50_et_on
>::digits
/ 2 - 2;
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
);
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
;
505 show_minima
<cpp_dec_float_50_et_on
>();
510 std::cout
.precision(std::numeric_limits
<cpp_dec_float_50_et_off
>::digits10
);
512 int bits
= std::numeric_limits
<cpp_dec_float_50_et_off
>::digits
/ 2 - 2;
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");
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
);
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
;
526 show_minima
<cpp_dec_float_50_et_off
>();
535 Typical output MSVC 15.7.3
537 brent_minimise_example.cpp
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.
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
557 Type double with limited iterations.
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.
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
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.
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.
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
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
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.
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
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
632 ============================================================================================================
634 // GCC 7.2.0 with quadmath
636 Brent's minimisation examples.
638 Type double - unlimited iterations (unwise?)
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
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
646 Type double with limited iterations.
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.
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
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.
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.
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.
671 x == 1 (compared to uncertainty 0.00034527) is true
672 f(x) == (0 compared to uncertainty 0.00034527) is true
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.
681 x == 1 (compared to uncertainty 1.490116e-008) is true
682 f(x) == (0 compared to uncertainty 1.490116e-008) is true
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.
691 x == 1 (compared to uncertainty 3.292723e-010) is true
692 f(x) == (0 compared to uncertainty 3.292723e-010) is true
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.
701 x == 1 (compared to uncertainty 1.38777878e-17) is true
702 f(x) == (0 compared to uncertainty 1.38777878e-17) is true
703 Bracketing -4.0000000000000000000000000000000000000000000000000 to 1.3333333333333333333333333333333333333333333333333
704 x at minimum = 0.99999999999999999999999999998813903221565569205253,
705 f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58, after 14 iterations.
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
717 x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58
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