> Can you think of a relatively safe way to get C++ to play nicely with
> the evaluation of these "tiny imaginary components are OK" problem? I
> get that a^b is "problematic" (to say the least) for a complex and b
> non-integer, but that's the problem I'm banging my head against at the
> moment.

Are you actually concerned about the non-integer case? In general, I don't think you're going to be able to use anything less than a true automatic differentiation approach there. The accuracy of the complex step derivative technique (which must of course be perfect to use imaginary components smaller than epsilon times the real components) depends absolutely on every function involved being exactly real or imaginary for a real input, and (the principal value of) the power function for _any_ non-integer exponent other than +/-1/2 doesn't satisfy that condition (mathematically, let alone in floating-point).

That said, for any _root_ that is 2 or odd, the problem is solvable. (In fact, std::sqrt(std::complex<...>) seems to already have the desired behavior in libstdc++ but not in libc++.) When the real component is negative, negate the entire number before and after applying std::pow, so that the conversion to polar form preserves the separation of the tiny imaginary component. This isn't quite a drop-in replacement for generic code, but it's probably a bit easier than real AD or your own custom "power" function that explicitly evaluates the imaginary component separately.

Davis

PS - Roots of course have the potential problem of having unbounded derivatives near 0, which this technique represents as the large value h^(1/p-1) sin(pi/2p).