Complex Step Differentiation

complex_step.cc

#include <stdio.h>
#include <complex>
#include <tuple>

namespace {

// See http://blogs.mathworks.com/cleve/2013/10/14/complex-step-differentiation/
template<class R, class F>
std::tuple<R, R> dF(F f, R x, R h=1e-8) {
  using C = std::complex<R>;
  C y = f(C(x, h));
  return std::make_tuple(y.real(), y.imag() / h);
}

}  // namespace

using namespace std;

using R = double;
using C = complex<R>;

C f(C x) {
  return exp(x) / (pow(cos(x), 3) + pow(sin(x), 3));
}

int main() {
  R x = M_PI / 4;
  R y, yd;
  tie(y, yd) = dF(f, x);
  printf("f(%f) = %f\tf'(%f) = %f\n", x, y, x, yd);
  return 0;
}