ceph/src/boost/libs/math/example/root_finding_multiprecision_example.cpp

   1 // Copyright Paul A. Bristow 2015.
   2
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 // Note that this file contains Quickbook mark-up as well as code
   9 // and comments, don't change any of the special comment mark-ups!
  10
  11 // Example of root finding using Boost.Multiprecision.
  12
  13 #ifndef BOOST_MATH_STANDALONE
  14
  15 #include <boost/math/tools/roots.hpp>
  16 //using boost::math::policies::policy;
  17 //using boost::math::tools::newton_raphson_iterate;
  18 //using boost::math::tools::halley_iterate;
  19 //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
  20 //using boost::math::tools::bracket_and_solve_root;
  21 //using boost::math::tools::toms748_solve;
  22
  23 #include <boost/math/special_functions/next.hpp> // For float_distance.
  24 #include <boost/math/special_functions/pow.hpp>
  25 #include <boost/math/constants/constants.hpp>
  26
  27 //[root_finding_multiprecision_include_1
  28 #include <boost/multiprecision/cpp_bin_float.hpp> // For cpp_bin_float_50.
  29 #include <boost/multiprecision/cpp_dec_float.hpp> // For cpp_dec_float_50.
  30 #ifndef _MSC_VER  // float128 is not yet supported by Microsoft compiler at 2013.
  31 #  include <boost/multiprecision/float128.hpp> // Requires libquadmath.
  32 #endif
  33 //] [/root_finding_multiprecision_include_1]
  34
  35 #include <iostream>
  36 // using std::cout; using std::endl;
  37 #include <iomanip>
  38 // using std::setw; using std::setprecision;
  39 #include <limits>
  40 // using std::numeric_limits;
  41 #include <tuple>
  42 #include <utility> // pair, make_pair
  43
  44 // #define BUILTIN_POW_GUESS // define to use std::pow function to obtain a guess.
  45
  46 template <class T>
  47 T cbrt_2deriv(T x)
  48 { // return cube root of x using 1st and 2nd derivatives and Halley.
  49   using namespace std;  // Help ADL of std functions.
  50   using namespace boost::math::tools; // For halley_iterate.
  51
  52   // If T is not a binary floating-point type, for example, cpp_dec_float_50
  53   // then frexp may not be defined,
  54   // so it may be necessary to compute the guess using a built-in type,
  55   // probably quickest using double, but perhaps with float or long double.
  56   // Note that the range of exponent may be restricted by a built-in-type for guess.
  57
  58   typedef long double guess_type;
  59
  60 #ifdef BUILTIN_POW_GUESS
  61   guess_type pow_guess = std::pow(static_cast<guess_type>(x), static_cast<guess_type>(1) / 3);
  62   T guess = pow_guess;
  63   T min = pow_guess /2;
  64   T max = pow_guess * 2;
  65 #else
  66   int exponent;
  67   frexp(static_cast<guess_type>(x), &exponent); // Get exponent of z (ignore mantissa).
  68   T guess = ldexp(static_cast<guess_type>(1.), exponent / 3); // Rough guess is to divide the exponent by three.
  69   T min = ldexp(static_cast<guess_type>(1.) / 2, exponent / 3); // Minimum possible value is half our guess.
  70   T max = ldexp(static_cast<guess_type>(2.), exponent / 3); // Maximum possible value is twice our guess.
  71 #endif
  72
  73   int digits = std::numeric_limits<T>::digits / 2; // Half maximum possible binary digits accuracy for type T.
  74   const std::uintmax_t maxit = 20;
  75   std::uintmax_t it = maxit;
  76   T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, it);
  77   // Can show how many iterations (updated by halley_iterate).
  78   // std::cout << "Iterations " << it << " (from max of "<< maxit << ")." << std::endl;
  79   return result;
  80 } // cbrt_2deriv(x)
  81
  82
  83 template <class T>
  84 struct cbrt_functor_2deriv
  85 { // Functor returning both 1st and 2nd derivatives.
  86   cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  87   { // Constructor stores value to find root of, for example:
  88   }
  89
  90   // using boost::math::tuple; // to return three values.
  91   std::tuple<T, T, T> operator()(T const& x)
  92   {
  93     // Return both f(x) and f'(x) and f''(x).
  94     T fx = x*x*x - a;                     // Difference (estimate x^3 - value).
  95     // std::cout << "x = " << x << "\nfx = " << fx << std::endl;
  96     T dx = 3 * x*x;                       // 1st derivative = 3x^2.
  97     T d2x = 6 * x;                        // 2nd derivative = 6x.
  98     return std::make_tuple(fx, dx, d2x);  // 'return' fx, dx and d2x.
  99   }
 100 private:
 101   T a;                                    // to be 'cube_rooted'.
 102 }; // struct cbrt_functor_2deriv
 103
 104 template <int n, class T>
 105 struct nth_functor_2deriv
 106 { // Functor returning both 1st and 2nd derivatives.
 107
 108   nth_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of)
 109   { /* Constructor stores value to find root of, for example: */ }
 110
 111   // using std::tuple; // to return three values.
 112   std::tuple<T, T, T> operator()(T const& x)
 113   {
 114     // Return both f(x) and f'(x) and f''(x).
 115     using boost::math::pow;
 116     T fx = pow<n>(x) - value;              // Difference (estimate x^3 - value).
 117     T dx = n * pow<n - 1>(x);              // 1st derivative = 5x^4.
 118     T d2x = n * (n - 1) * pow<n - 2 >(x);  // 2nd derivative = 20 x^3
 119     return std::make_tuple(fx, dx, d2x);   // 'return' fx, dx and d2x.
 120   }
 121 private:
 122   T value;                                 // to be 'nth_rooted'.
 123 }; // struct nth_functor_2deriv
 124
 125
 126 template <int n, class T>
 127 T nth_2deriv(T x)
 128 {
 129   // return nth root of x using 1st and 2nd derivatives and Halley.
 130   using namespace std;  // Help ADL of std functions.
 131   using namespace boost::math; // For halley_iterate.
 132
 133   int exponent;
 134   frexp(x, &exponent);                                 // Get exponent of z (ignore mantissa).
 135   T guess = ldexp(static_cast<T>(1.), exponent / n);   // Rough guess is to divide the exponent by three.
 136   T min = ldexp(static_cast<T>(0.5), exponent / n);    // Minimum possible value is half our guess.
 137   T max = ldexp(static_cast<T>(2.), exponent / n);     // Maximum possible value is twice our guess.
 138
 139   int digits = std::numeric_limits<T>::digits / 2;     // Half maximum possible binary digits accuracy for type T.
 140   const std::uintmax_t maxit = 50;
 141   std::uintmax_t it = maxit;
 142   T result = halley_iterate(nth_functor_2deriv<n, T>(x), guess, min, max, digits, it);
 143   // Can show how many iterations (updated by halley_iterate).
 144   std::cout << it << " iterations (from max of " << maxit << ")" << std::endl;
 145
 146   return result;
 147 } // nth_2deriv(x)
 148
 149 //[root_finding_multiprecision_show_1
 150
 151 template <typename T>
 152 T show_cube_root(T value)
 153 { // Demonstrate by printing the root using all definitely significant digits.
 154   std::cout.precision(std::numeric_limits<T>::digits10);
 155   T r = cbrt_2deriv(value);
 156   std::cout << "value = " << value << ", cube root =" << r << std::endl;
 157   return r;
 158 }
 159
 160 //] [/root_finding_multiprecision_show_1]
 161
 162 int main()
 163 {
 164   std::cout << "Multiprecision Root finding Example." << std::endl;
 165   // Show all possibly significant decimal digits.
 166   std::cout.precision(std::numeric_limits<double>::digits10);
 167   // or use   cout.precision(max_digits10 = 2 + std::numeric_limits<double>::digits * 3010/10000);
 168   //[root_finding_multiprecision_example_1
 169   using boost::multiprecision::cpp_dec_float_50; // decimal.
 170   using boost::multiprecision::cpp_bin_float_50; // binary.
 171 #ifndef _MSC_VER  // Not supported by Microsoft compiler.
 172   using boost::multiprecision::float128;
 173 #endif
 174   //] [/root_finding_multiprecision_example_1
 175
 176   try
 177   { // Always use try'n'catch blocks with Boost.Math to get any error messages.
 178     // Increase the precision to 50 decimal digits using Boost.Multiprecision
 179 //[root_finding_multiprecision_example_2
 180
 181       std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
 182
 183       cpp_dec_float_50 two = 2; //
 184       cpp_dec_float_50  r = cbrt_2deriv(two);
 185       std::cout << "cbrt(" << two << ") = " << r << std::endl;
 186
 187       r = cbrt_2deriv(2.); // Passing a double, so ADL will compute a double precision result.
 188       std::cout << "cbrt(" << two << ") = " << r << std::endl;
 189       // cbrt(2) = 1.2599210498948731906665443602832965552806854248047 'wrong' from digits 17 onwards!
 190       r = cbrt_2deriv(static_cast<cpp_dec_float_50>(2.)); // Passing a cpp_dec_float_50,
 191       // so will compute a cpp_dec_float_50 precision result.
 192       std::cout << "cbrt(" << two << ") = " << r << std::endl;
 193       r = cbrt_2deriv<cpp_dec_float_50>(2.); // Explicitly a cpp_dec_float_50, so will compute a cpp_dec_float_50 precision result.
 194       std::cout << "cbrt(" << two << ") = " << r << std::endl;
 195       // cpp_dec_float_50 1.2599210498948731647672106072782283505702514647015
 196 //] [/root_finding_multiprecision_example_2
 197      //  N[2^(1/3), 50]  1.2599210498948731647672106072782283505702514647015
 198
 199       //show_cube_root(2); // Integer parameter - Errors!
 200       //show_cube_root(2.F); // Float parameter - Warnings!
 201 //[root_finding_multiprecision_example_3
 202       show_cube_root(2.);
 203       show_cube_root(2.L);
 204       show_cube_root(two);
 205
 206 //] [/root_finding_multiprecision_example_3
 207
 208   }
 209   catch (const std::exception& e)
 210   { // Always useful to include try&catch blocks because default policies
 211     // are to throw exceptions on arguments that cause errors like underflow & overflow.
 212     // Lacking try&catch blocks, the program will abort without a message below,
 213     // which may give some helpful clues as to the cause of the exception.
 214     std::cout <<
 215       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
 216   }
 217   return 0;
 218 } // int main()
 219
 220
 221 /*
 222
 223 Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_multiprecision.exe"
 224 Multiprecision Root finding Example.
 225 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
 226 cbrt(2) = 1.2599210498948731906665443602832965552806854248047
 227 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
 228 cbrt(2) = 1.2599210498948731647672106072782283505702514647015
 229 value = 2, cube root =1.25992104989487
 230 value = 2, cube root =1.25992104989487
 231 value = 2, cube root =1.2599210498948731647672106072782283505702514647015
 232
 233
 234 */
 235
 236 #endif // BOOST_MATH_STANDALONE