Dual will not work correctly for complex numbers

[petar-andrejic]

As it stands the following overloads seem to be problematic, and restrict the type to work only for real types:

template<typename R, Requires<isExpr<R>> = true> constexpr auto abs2(R&& r) { return std::forward<R>(r) * std::forward<R>(r); }
template<typename R, Requires<isExpr<R>> = true> constexpr auto conj(R&& r) { return std::forward<R>(r); }
template<typename R, Requires<isExpr<R>> = true> constexpr auto real(R&& r) { return std::forward<R>(r); }
template<typename R, Requires<isExpr<R>> = true> constexpr auto imag(R&&) { return 0.0; }

This contradicts the info given in the tutorial, which suggests using Dual<std::complex<double>,std::complex<double>>

A simple example of how it goes wrong:

#include <complex>
#include <iostream>

#include <autodiff/forward/dual.hpp>

using namespace autodiff;
using namespace std::complex_literals;
using cxdouble = std::complex<double>;
using cxdual = Dual<cxdouble, cxdouble>;

namespace autodiff::detail
{

template <typename T> struct ArithmeticTraits<std::complex<T>> : ArithmeticTraits<T>
{
};
} 

int main()
{
    cxdual x = 1.0 + 0.3i;
    std::cout << "x = 1.0 + 0.3i \n";
    std::cout << "real(x) = " << real(x).val << "\n"; // this should be a double, not a complex
    std::cout << "imag(x) = " << imag(x).val << "\n"; // this shouldn't be zero but it is!
    std::cout << "conj(x) = " << conj(x).val.real() << " + " << conj(x).val.imag() << "i\n"; // this shouldn't be the same as x!
}

The output should be ideally

x = 1.0 + 0.3i
real(x) = 1
imag(x) = 0.3
conj(x) = 1 + -0.3i

or similar
But instead we get

x = 1.0 + 0.3i 
real(x) = (1,0.3)
imag(x) = 0
conj(x) = 1 + 0.3i

which matches what we can see the hardcoded methods are actually doing. Abs2 is doing a similarly problematic thing by returning the number times itself no matter what type we are passing in.

As it stands, it should be fine to use std::real, std::imag, not so sure about std::conj. The first two work correctly for real types anyway, while std::conj on a double returns a complex which could be a problem.

For abs2 and abs, this would have to be changed as well, abs2 to work correctly even for complex types, and in the apply method for abs binary comparisons are used, which will not work for complex numbers

[petar-andrejic]

I'll try have a look at implementing the changes myself, and will open a pull request if I succeed

[allanleal]

Hi @petar-andrejic , thanks for initiating this discussion. I don't deal with complex numbers and I'm currently very busy with my main project (Reaktoro - https://reaktoro.org) which uses autodiff::real. Thus, if you could find a solution for complex numbers with dual and provide a Pull Request, that would be great. Thanks!

[petar-andrejic]

So digging a bit further, the main issue is that the current implementation only makes sense for holomorphic complex functions, since otherwise one has to consider wirtinger derivatives (which would have to be a new api specifically for complex numbers and as such probably out of scope since you're busy with other things). The current implementation works fine as long as one uses only holomorphic functions, but things like real, imag, abs and conj are not holomorphic, and ideally they would somehow not be implemented and refuse to compile for Dual<complex, complex>, I'm not sure what a clean way of doing that would be though.

As it currently stands the implementation should still be fine for std::complex, i.e. using std::complex to take the derivative of complex valued functions of a real variable. However, that doesn't work for me since the compiler doesn't find an appropriate template, e.g.

#include <complex>
#include <iostream>

#include <autodiff/forward/dual.hpp>

using namespace autodiff;
using namespace std::complex_literals;
using cxdual = std::complex<dual>;

const cxdual im = cxdual(1i);

int main()
{
    dual w = 1.0;
    cxdual phi = im * w;
    cxdual f = exp(phi);
    std::cout << f.real().val << "\n";
    return 0;
}

the above won't compile since it can't find an appropriate template for cos and sin