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 m_n 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(dx⁶), 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(dx⁶). 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;