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

   1 // Copyright Christopher Kormanyos 2013.
   2 // Distributed under the Boost Software License, Version 1.0.
   3 // (See accompanying file LICENSE_1_0.txt or
   4 // copy at http://www.boost.org/LICENSE_1_0.txt).
   5
   6 #ifdef _MSC_VER
   7 #  pragma warning (disable : 4996) // assignment operator could not be generated.
   8 #endif
   9
  10 # include <iostream>
  11 # include <iomanip>
  12 # include <limits>
  13 # include <cmath>
  14
  15 #include <boost/static_assert.hpp>
  16 #include <boost/type_traits/is_floating_point.hpp>
  17 #include <boost/math/special_functions/next.hpp> // for float_distance
  18
  19 //[numeric_derivative_example
  20 /*`The following example shows how multiprecision calculations can be used to
  21 obtain full precision in a numerical derivative calculation that suffers from precision loss.
  22
  23 Consider some well-known central difference rules for numerically
  24 computing the 1st derivative of a function [f'(x)] with [/x] real.
  25
  26 Need a reference here?  Introduction to Partial Differential Equations, Peter J. Olver
  27  December 16, 2012
  28
  29 Here, the implementation uses a C++ template that can be instantiated with various
  30 floating-point types such as `float`, `double`, `long double`, or even
  31 a user-defined floating-point type like __multiprecision.
  32
  33 We will now use the derivative template with the built-in type `double` in
  34 order to numerically compute the derivative of a function, and then repeat
  35 with a 5 decimal digit higher precision user-defined floating-point type.
  36
  37 Consider the function  shown below.
  38 !!
  39 (3)
  40 We will now take the derivative of this function with respect to x evaluated
  41 at x = 3= 2. In other words,
  42
  43 (4)
  44
  45 The expected result is
  46
  47  0:74535 59924 99929 89880 . (5)
  48 The program below uses the derivative template in order to perform
  49 the numerical calculation of this derivative. The program also compares the
  50 numerically-obtained result with the expected result and reports the absolute
  51 relative error scaled to a deviation that can easily be related to the number of
  52 bits of lost precision.
  53
  54 */
  55
  56 /*` [note Requires the C++11 feature of
  57 [@http://en.wikipedia.org/wiki/Anonymous_function#C.2B.2B anonymous functions]
  58 for the derivative function calls like `[]( const double & x_) -> double`.
  59 */
  60
  61
  62
  63 template <typename value_type,  typename function_type>
  64 value_type derivative (const value_type x, const value_type dx, function_type function)
  65 {
  66   /*! \brief Compute the derivative of function using a 3-point central difference rule of O(dx^6).
  67     \tparam value_type, floating-point type, for example: `double` or `cpp_dec_float_50`
  68     \tparam function_type
  69
  70     \param x Value at which to evaluate derivative.
  71     \param dx Incremental step-size.
  72     \param function Function whose derivative is to computed.
  73
  74     \return derivative at x.
  75   */
  76
  77   BOOST_STATIC_ASSERT_MSG(false == std::numeric_limits<value_type>::is_integer, "value_type must be a floating-point type!");
  78
  79   const value_type dx2(dx * 2U);
  80   const value_type dx3(dx * 3U);
  81   // Difference terms.
  82   const value_type m1 ((function (x + dx) - function(x - dx)) / 2U);
  83   const value_type m2 ((function (x + dx2) - function(x - dx2)) / 4U);
  84   const value_type m3 ((function (x + dx3) - function(x - dx3)) / 6U);
  85   const value_type fifteen_m1 (m1 * 15U);
  86   const value_type six_m2 (m2 * 6U);
  87   const value_type ten_dx (dx * 10U);
  88   return ((fifteen_m1 - six_m2) + m3) / ten_dx;  // Derivative.
  89 } //
  90
  91 #include <boost/multiprecision/cpp_dec_float.hpp>
  92   using boost::multiprecision::number;
  93   using boost::multiprecision::cpp_dec_float;
  94
  95 // Re-compute using 5 extra decimal digits precision (22) than double (17).
  96 #define MP_DIGITS10 unsigned (std::numeric_limits<double>::max_digits10 + 5)
  97
  98 typedef cpp_dec_float<MP_DIGITS10> mp_backend;
  99 typedef number<mp_backend> mp_type;
 100
 101
 102 int main()
 103 {
 104   {
 105     const double d =
 106       derivative
 107       ( 1.5, // x = 3.2
 108         std::ldexp (1., -9), // step size 2^-9 = see below for choice.
 109         [](const double & x)->double // Function f(x).
 110         {
 111           return std::sqrt((x * x) - 1.) - std::acos(1. / x);
 112         }
 113       );
 114
 115     // The 'exactly right' result is [sqrt]5 / 3 = 0.74535599249992989880.
 116     const double rel_error = (d - 0.74535599249992989880) / 0.74535599249992989880;
 117     const double bit_error = std::abs(rel_error) / std::numeric_limits<double>::epsilon();
 118     std::cout.precision (std::numeric_limits<double>::digits10); // Show all guaranteed decimal digits.
 119     std::cout << std::showpoint ; // Ensure that any trailing zeros are shown too.
 120
 121     std::cout << " derivative : " << d << std::endl;
 122     std::cout << " expected   : " << 0.74535599249992989880 << std::endl;
 123     // Can compute an 'exact' value using multiprecision type.
 124     std::cout << " expected   : " << sqrt(static_cast<mp_type>(5))/3U << std::endl;
 125     std::cout << " bit_error : " << static_cast<unsigned long>(bit_error)  << std::endl;
 126
 127     std::cout.precision(6);
 128     std::cout << "float_distance = " << boost::math::float_distance(0.74535599249992989880, d) << std::endl;
 129
 130   }
 131
 132   { // Compute using multiprecision type with an extra 5 decimal digits of precision.
 133     const mp_type mp =
 134       derivative(mp_type(mp_type(3) / 2U), // x = 3/2
 135         mp_type(mp_type(1) / 10000000U), // Step size 10^7.
 136         [](const mp_type & x)->mp_type
 137         {
 138           return sqrt((x * x) - 1.) - acos (1. / x); // Function
 139         }
 140     );
 141
 142     const double d = mp.convert_to<double>(); // Convert to closest double.
 143     const double rel_error = (d - 0.74535599249992989880) / 0.74535599249992989880;
 144     const double bit_error = std::abs (rel_error) / std::numeric_limits<double>::epsilon();
 145     std::cout.precision (std::numeric_limits <double>::digits10); // All guaranteed decimal digits.
 146     std::cout << std::showpoint ; // Ensure that any trailing zeros are shown too.
 147     std::cout << " derivative : " << d << std::endl;
 148     // Can compute an 'exact' value using multiprecision type.
 149     std::cout << " expected   : " << sqrt(static_cast<mp_type>(5))/3U << std::endl;
 150     std::cout << " expected   : " << 0.74535599249992989880
 151     << std::endl;
 152     std::cout << " bit_error : "  << static_cast<unsigned long>(bit_error)  << std::endl;
 153
 154     std::cout.precision(6);
 155     std::cout << "float_distance = " << boost::math::float_distance(0.74535599249992989880, d) << std::endl;
 156
 157
 158   }
 159
 160
 161 } // int main()
 162
 163 /*`
 164 The result of this program on a system with an eight-byte, 64-bit IEEE-754
 165 conforming floating-point representation for `double` is:
 166
 167  derivative : 0.745355992499951
 168
 169  derivative : 0.745355992499943
 170  expected   : 0.74535599249993
 171  bit_error : 78
 172
 173     derivative : 0.745355992499930
 174    expected   : 0.745355992499930
 175    bit_error : 0
 176
 177 The resulting bit error is 0. This means that the result of the derivative
 178 calculation is bit-identical with the double representation of the expected result,
 179 and this is the best result possible for the built-in type.
 180
 181 The derivative in this example has a known closed form. There are, however,
 182 countless situations in numerical analysis (and not only for numerical deriva-
 183 tives) for which the calculation at hand does not have a known closed-form
 184 solution or for which the closed-form solution is highly inconvenient to use. In
 185 such cases, this technique may be useful.
 186
 187 This example has shown how multiprecision can be used to add extra digits
 188 to an ill-conditioned calculation that suffers from precision loss. When the result
 189 of the multiprecision calculation is converted to a built-in type such as double,
 190 the entire precision of the result in double is preserved.
 191
 192  */
 193
 194 /*
 195
 196   Description: Autorun "J:\Cpp\big_number\Debug\numerical_derivative_example.exe"
 197    derivative : 0.745355992499943
 198    expected   : 0.745355992499930
 199    expected   : 0.745355992499930
 200    bit_error : 78
 201   float_distance = 117.000
 202    derivative : 0.745355992499930
 203    expected   : 0.745355992499930
 204    expected   : 0.745355992499930
 205    bit_error : 0
 206   float_distance = 0.000000
 207
 208  */
 209