The recent "rediscovery" of complex step derivatives (https://sinews.siam.org/Details-Page/differentiation-without-a-difference) has made numerical differentiation as accurate as any other evaluation in double precision arithmetic. To fully make use of this technique, all functions must accept either complex or double arguments. In principle that is no problem for C++. In practice, serious problems occur in some cases.

Here first is a simple example in Python of when things go right. The derivative of x^2.0 is 2.0*x, so the derivative of x^2.0 at x=-0.1 should be dy/dx=-0.2. In Python, no problem to use complex step derivatives to evaluate:


h = 1e-100

z = -0.1+1j*h

print(pow(z, 2.0), pow(z, 2.0).imag/h, (z*z).imag/h)

gives


(0.010000000000000002-2e-101j) -0.2 -0.2

On the contrary, the same example in C++


#include <iostream>

#include <complex>

#include <cmath>



int main()

{

    std::cout << pow(std::complex<double>(-0.1, 1e-100), 2.0) << std::endl;

}

gives


(0.01,-2.44929e-18)


I believe the problem has to do with the handling of the branch-cut of the log function. In any case, this demonstrates a result that is silently in error by 83 orders of magnitude! Had I multiplied the complex step by itself rather than pow(z,2.0), I would have obtained the correct result.

I realize that I am probing an uncomfortable part of the complex plane, but I wonder if this could be handled more like Python, to minimize surprises for complex step derivative approaches?

Ian