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

   1 // Copyright Paul A. Bristow 2014, 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 finding nth root using 1st and 2nd derivatives of x^n.
  12
  13 #include <boost/math/tools/roots.hpp>
  14 //using boost::math::policies::policy;
  15 //using boost::math::tools::newton_raphson_iterate;
  16 //using boost::math::tools::halley_iterate;
  17 //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
  18 //using boost::math::tools::bracket_and_solve_root;
  19 //using boost::math::tools::toms748_solve;
  20
  21 #include <boost/math/special_functions/next.hpp>
  22 #include <boost/multiprecision/cpp_dec_float.hpp>
  23 #include <boost/math/special_functions/pow.hpp>
  24 #include <boost/math/constants/constants.hpp>
  25
  26 #include <boost/multiprecision/cpp_dec_float.hpp> // For cpp_dec_float_50.
  27 #include <boost/multiprecision/cpp_bin_float.hpp> // using boost::multiprecision::cpp_bin_float_50;
  28 #ifndef _MSC_VER  // float128 is not yet supported by Microsoft compiler at 2013.
  29 # include <boost/multiprecision/float128.hpp>
  30 #endif
  31
  32 #include <iostream>
  33 // using std::cout; using std::endl;
  34 #include <iomanip>
  35 // using std::setw; using std::setprecision;
  36 #include <limits>
  37 using std::numeric_limits;
  38 #include <tuple>
  39 #include <utility> // pair, make_pair
  40
  41
  42 //[root_finding_nth_functor_2deriv
  43 template <int N, class T = double>
  44 struct nth_functor_2deriv
  45 { // Functor returning both 1st and 2nd derivatives.
  46   static_assert(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
  47   static_assert((N > 0) == true, "root N must be > 0!");
  48
  49   nth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  50   { /* Constructor stores value a to find root of, for example: */ }
  51
  52   // using boost::math::tuple; // to return three values.
  53   std::tuple<T, T, T> operator()(T const& x)
  54   {
  55     // Return f(x), f'(x) and f''(x).
  56     using boost::math::pow;
  57     T fx = pow<N>(x) - a;                  // Difference (estimate x^n - a).
  58     T dx = N * pow<N - 1>(x);              // 1st derivative f'(x).
  59     T d2x = N * (N - 1) * pow<N - 2 >(x);  // 2nd derivative f''(x).
  60
  61     return std::make_tuple(fx, dx, d2x);   // 'return' fx, dx and d2x.
  62   }
  63 private:
  64   T a;                                     // to be 'nth_rooted'.
  65 };
  66
  67 //] [/root_finding_nth_functor_2deriv]
  68
  69 /*
  70 To show the progress, one might use this before the return statement above?
  71 #ifdef BOOST_MATH_ROOT_DIAGNOSTIC
  72 std::cout << " x = " << x << ", fx = " << fx << ", dx = " << dx << ", dx2 = " << d2x << std::endl;
  73 #endif
  74 */
  75
  76 // If T is a floating-point type, might be quicker to compute the guess using a built-in type,
  77 // probably quickest using double, but perhaps with float or long double, T.
  78
  79 // If T is a type for which frexp and ldexp are not defined,
  80 // then it is necessary to compute the guess using a built-in type,
  81 // probably quickest (but limited range) using double,
  82 // but perhaps with float or long double, or a multiprecision T for the full range of T.
  83 // typedef double guess_type; is used to specify the this.
  84
  85 //[root_finding_nth_function_2deriv
  86
  87 template <int N, class T = double>
  88 T nth_2deriv(T x)
  89 { // return nth root of x using 1st and 2nd derivatives and Halley.
  90
  91   using namespace std;  // Help ADL of std functions.
  92   using namespace boost::math::tools; // For halley_iterate.
  93
  94   static_assert(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
  95   static_assert((N > 0) == true, "root N must be > 0!");
  96   static_assert((N > 1000) == false, "root N is too big!");
  97
  98   typedef double guess_type; // double may restrict (exponent) range for a multiprecision T?
  99
 100   int exponent;
 101   frexp(static_cast<guess_type>(x), &exponent);                 // Get exponent of z (ignore mantissa).
 102   T guess = ldexp(static_cast<guess_type>(1.), exponent / N);   // Rough guess is to divide the exponent by n.
 103   T min = ldexp(static_cast<guess_type>(1.) / 2, exponent / N); // Minimum possible value is half our guess.
 104   T max = ldexp(static_cast<guess_type>(2.), exponent / N);     // Maximum possible value is twice our guess.
 105
 106   int digits = std::numeric_limits<T>::digits * 0.4;            // Accuracy triples with each step, so stop when
 107                                                                 // slightly more than one third of the digits are correct.
 108   const std::uintmax_t maxit = 20;
 109   std::uintmax_t it = maxit;
 110   T result = halley_iterate(nth_functor_2deriv<N, T>(x), guess, min, max, digits, it);
 111   return result;
 112 }
 113
 114 //] [/root_finding_nth_function_2deriv]
 115
 116
 117 template <int N, typename T = double>
 118 T show_nth_root(T value)
 119 { // Demonstrate by printing the nth root using all possibly significant digits.
 120   //std::cout.precision(std::numeric_limits<T>::max_digits10);
 121   // or use   cout.precision(max_digits10 = 2 + std::numeric_limits<double>::digits * 3010/10000);
 122   // Or guaranteed significant digits:
 123    std::cout.precision(std::numeric_limits<T>::digits10);
 124
 125   T r = nth_2deriv<N>(value);
 126   std::cout << "Type " << typeid(T).name() << " value = " << value << ", " << N << "th root = " << r << std::endl;
 127   return r;
 128 } // print_nth_root
 129
 130
 131 int main()
 132 {
 133   std::cout << "nth Root finding Example." << std::endl;
 134   using boost::multiprecision::cpp_dec_float_50; // decimal.
 135   using boost::multiprecision::cpp_bin_float_50; // binary.
 136 #ifndef _MSC_VER  // Not supported by Microsoft compiler.
 137   using boost::multiprecision::float128; // Requires libquadmath
 138 #endif
 139   try
 140   { // Always use try'n'catch blocks with Boost.Math to get any error messages.
 141
 142 //[root_finding_n_example_1
 143     double r1 = nth_2deriv<5, double>(2); // Integral value converted to double.
 144
 145     // double r2 = nth_2deriv<5>(2); // Only floating-point type types can be used!
 146
 147 //] [/root_finding_n_example_1
 148
 149     //show_nth_root<5, float>(2); // Integral value converted to float.
 150     //show_nth_root<5, float>(2.F); // 'initializing' : conversion from 'double' to 'float', possible loss of data
 151
 152 //[root_finding_n_example_2
 153
 154
 155     show_nth_root<5, double>(2.);
 156     show_nth_root<5, long double>(2.);
 157 #ifndef _MSC_VER  // float128 is not supported by Microsoft compiler 2013.
 158     show_nth_root<5, float128>(2);
 159 #endif
 160     show_nth_root<5, cpp_dec_float_50>(2); // dec
 161     show_nth_root<5, cpp_bin_float_50>(2); // bin
 162 //] [/root_finding_n_example_2
 163
 164     // show_nth_root<1000000>(2.); // Type double value = 2, 555th root = 1.00124969405651
 165     // Type double value = 2, 1000th root = 1.00069338746258
 166     // Type double value = 2, 1000000th root = 1.00000069314783
 167   }
 168   catch (const std::exception& e)
 169   { // Always useful to include try & catch blocks because default policies
 170     // are to throw exceptions on arguments that cause errors like underflow, overflow.
 171     // Lacking try & catch blocks, the program will abort without a message below,
 172     // which may give some helpful clues as to the cause of the exception.
 173     std::cout <<
 174       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
 175   }
 176   return 0;
 177 } // int main()
 178
 179
 180 /*
 181 //[root_finding_example_output_1
 182  Using MSVC 2013
 183
 184 nth Root finding Example.
 185 Type double value = 2, 5th root = 1.14869835499704
 186 Type long double value = 2, 5th root = 1.14869835499704
 187 Type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_dec_float<50,int,void>,1> value = 2,
 188   5th root = 1.1486983549970350067986269467779275894438508890978
 189 Type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,0> value = 2,
 190   5th root = 1.1486983549970350067986269467779275894438508890978
 191
 192 //] [/root_finding_example_output_1]
 193
 194 //[root_finding_example_output_2
 195
 196  Using GCC 4.91  (includes float_128 type)
 197
 198  nth Root finding Example.
 199 Type d value = 2, 5th root = 1.14869835499704
 200 Type e value = 2, 5th root = 1.14869835499703501
 201 Type N5boost14multiprecision6numberINS0_8backends16float128_backendELNS0_26expression_template_optionE0EEE value = 2, 5th root = 1.148698354997035006798626946777928
 202 Type N5boost14multiprecision6numberINS0_8backends13cpp_dec_floatILj50EivEELNS0_26expression_template_optionE1EEE value = 2, 5th root = 1.1486983549970350067986269467779275894438508890978
 203 Type N5boost14multiprecision6numberINS0_8backends13cpp_bin_floatILj50ELNS2_15digit_base_typeE10EviLi0ELi0EEELNS0_26expression_template_optionE0EEE value = 2, 5th root = 1.1486983549970350067986269467779275894438508890978
 204
 205 RUN SUCCESSFUL (total time: 63ms)
 206
 207 //] [/root_finding_example_output_2]
 208 */
 209
 210 /*
 211 Throw out of range using GCC release mode :-(
 212
 213  */