...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
In this example we'll add even more power to generic numeric programming using not only different floating-point types but also function objects as template parameters. Consider some well-known central difference rules for numerically computing the first derivative of a function f′(x) with x ∈ ℜ:
Where the difference terms mn are given by:
and dx is the step-size of the derivative.
The third formula in Equation 1 is a three-point central difference rule. It calculates the first derivative of f′(x) to O(dx6), where dx is the given step-size. For example, if the step-size is 0.01 this derivative calculation has about 6 decimal digits of precision - just about right for the 7 decimal digits of single-precision float. Let's make a generic template subroutine using this three-point central difference rule. In particular:
template<typename value_type, typename function_type> value_type derivative(const value_type x, const value_type dx, function_type func) { // Compute d/dx[func(*first)] using a three-point // central difference rule of O(dx^6). const value_type dx1 = dx; const value_type dx2 = dx1 * 2; const value_type dx3 = dx1 * 3; const value_type m1 = (func(x + dx1) - func(x - dx1)) / 2; const value_type m2 = (func(x + dx2) - func(x - dx2)) / 4; const value_type m3 = (func(x + dx3) - func(x - dx3)) / 6; const value_type fifteen_m1 = 15 * m1; const value_type six_m2 = 6 * m2; const value_type ten_dx1 = 10 * dx1; return ((fifteen_m1 - six_m2) + m3) / ten_dx1; }
The derivative()
template function can be used to compute the first derivative of any
function to O(dx6). For example, consider the first
derivative of sin(x) evaluated at x =
π/3. In other words,
The code below computes the derivative in Equation 3 for float, double and boost's multiple-precision type cpp_dec_float_50.
#include <iostream> #include <iomanip> #include <boost/multiprecision/cpp_dec_float.hpp> #include <boost/math/constants/constants.hpp> int main(int, char**) { using boost::math::constants::pi; using boost::multiprecision::cpp_dec_float_50; // // We'll pass a function pointer for the function object passed to derivative, // the typecast is needed to select the correct overload of std::sin: // const float d_f = derivative( pi<float>() / 3, 0.01F, static_cast<float(*)(float)>(std::sin) ); const double d_d = derivative( pi<double>() / 3, 0.001, static_cast<double(*)(double)>(std::sin) ); // // In the cpp_dec_float_50 case, the sin function is multiply overloaded // to handle expression templates etc. As a result it's hard to take its // address without knowing about its implementation details. We'll use a // C++11 lambda expression to capture the call. // We also need a typecast on the first argument so we don't accidentally pass // an expression template to a template function: // const cpp_dec_float_50 d_mp = derivative( cpp_dec_float_50(pi<cpp_dec_float_50>() / 3), cpp_dec_float_50(1.0E-9), [](const cpp_dec_float_50& x) -> cpp_dec_float_50 { return sin(x); } ); // 5.000029e-001 std::cout << std::setprecision(std::numeric_limits<float>::digits10) << d_f << std::endl; // 4.999999999998876e-001 std::cout << std::setprecision(std::numeric_limits<double>::digits10) << d_d << std::endl; // 4.99999999999999999999999999999999999999999999999999e-01 std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << d_mp << std::endl; }
The expected value of the derivative is 0.5. This central difference rule in this example is ill-conditioned, meaning it suffers from slight loss of precision. With that in mind, the results agree with the expected value of 0.5.
We can take this a step further and use our derivative function to compute a partial derivative. For example if we take the incomplete gamma function P(a, z), and take the derivative with respect to z at (2,2) then we can calculate the result as shown below, for good measure we'll compare with the "correct" result obtained from a call to gamma_p_derivative, the results agree to approximately 44 digits:
cpp_dec_float_50 gd = derivative( cpp_dec_float_50(2), cpp_dec_float_50(1.0E-9), [](const cpp_dec_float_50& x) ->cpp_dec_float_50 { return boost::math::gamma_p(2, x); } ); // 2.70670566473225383787998989944968806815263091819151e-01 std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << gd << std::endl; // 2.70670566473225383787998989944968806815253190143120e-01 std::cout << boost::math::gamma_p_derivative(cpp_dec_float_50(2), cpp_dec_float_50(2)) << std::endl;