libs/math/example/brent_minimise_example.cpp

   1 //! \file
   2 //! \brief Brent_minimise_example.cpp
   3
   4 // Copyright Paul A. Bristow 2015, 2018.
   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>
  28 #include <boost/test/tools/floating_point_comparison.hpp> // For is_close_at)tolerance and is_small
  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
  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>
  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
  61 // An example of a function for which we want to find a minimum.
  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)
  71   {
  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)
  82   {
  83     return (x + 3) * (x - 1) * (x - 1); // (x + 3)(x - 1)^2
  84   }
  85 };
  86 //] [/brent_minimise_T_functor]
  87
  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
  96 //[brent_minimise_close
  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)
 105 {
 106   using boost::math::fpc::close_at_tolerance;
 107   using boost::math::fpc::is_small;
 108   using boost::math::fpc::FPC_STRONG;
 109
 110   if (is_small<T>(expect, tolerance))
 111   {
 112     return is_small<T>(got, tolerance);
 113   }
 114
 115   return close_at_tolerance<T>(tolerance, FPC_STRONG) (expect, got);
 116 } // bool is_close(T expect, T got, T tolerance)
 117
 118 //] [/brent_minimise_close]
 119
 120 //[brent_minimise_T_show
 121
 122 //! Example template function to find and show minima.
 123 //! \tparam T floating-point or fixed_point type.
 124 template <class T>
 125 void show_minima()
 126 {
 127   using boost::math::tools::brent_find_minima;
 128   using std::sqrt;
 129   try
 130   { // Always use try'n'catch blocks with Boost.Math to ensure you get any error messages.
 131
 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.
 135
 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."
 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
 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);
 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());
 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.
 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 {
 189   using boost::math::tools::brent_find_minima;
 190   using std::sqrt;
 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.
 194
 195   // Tip - using
 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.
 200
 201   // Specific type double - unlimited iterations (unwise?).
 202   {
 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);
 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;
 224   using boost::math::fpc::is_small;
 225
 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]
 232
 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
 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())
 253     << std::endl;
 254
 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
 262   }
 263
 264   std::cout << "\nType double with limited iterations and half double bits." << std::endl;
 265   {
 266
 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.
 271
 272   std::cout << "Showing " << bits_div_2 << " bits precision with " << prec
 273     << " decimal digits from tolerance " << sqrt(epsilon_2)
 274     << std::endl;
 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.
 282
 283 //] [/brent_minimise_double_4]
 284   }
 285   // x at minimum = 1, f(1) = 5.04852568e-018
 286
 287   {
 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)
 295       << std::endl;
 296     std::streamsize precision_5 = std::cout.precision(prec); // Save & set.
 297     const boost::uintmax_t maxit = 20;
 298
 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.
 304
 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   {
 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;
 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;
 325     std::cout.precision(precision_t1);  // Restore.
 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>();
 336
 337  //] [/brent_minimise_template_fd]
 338
 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;
 342 #endif
 343 //] [/brent_minimise_mp_include_1]
 344
 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.
 348 #endif
 349 //] [/brent_minimise_template_quad
 350
 351   // User-defined floating-point template.
 352
 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.
 356
 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.
 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
 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);
 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");
 381
 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);
 387
 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;
 391
 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()));
 394
 395 //] [/brent_minimise_mp_1]
 396
 397 /*
 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
 411 14 iterations
 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
 448   // Some examples of switching expression templates on and off follow.
 449
 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
 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
 621
 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
 630
 631
 632 ============================================================================================================
 633
 634  // GCC 7.2.0 with quadmath
 635
 636 Brent's minimisation examples.
 637
 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
 645
 646 Type double with limited iterations.
 647 Precision bits = 53
 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
 655 10 iterations.
 656
 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.
 663
 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.
 671 x == 1 (compared to uncertainty 0.00034527) is true
 672 f(x) == (0 compared to uncertainty 0.00034527) is true
 673
 674
 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.
 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
 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.
 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
 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.
 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
 717 x at minimum = 0.99999999999999999999999999998813903221565569205253, f(0.99999999999999999999999999998813903221565569205253) = 5.6273022712501408640665300316078046703496236636624e-58
 718 14 iterations.
 719
 720
 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
 729
 730 */