Division and absolute value of std::complex<mpfr::mpreal> with undesired behavior

[weslleyspereira]

Hi. I am working on the <T>LAPACK project (https://github.com/tlapack/tlapack), which is an attempt to move the reference LAPACK to C++. We are trying to use the mpfr::mpreal datatype to enable any precision through this linear algebra library.
Related and unrelated, we were checking some behavior of various languages/compilers (C++, Fortran, etc.) for exception handling, and we are including MPFR in our tests as much as possible.

As a result of the tests with complex numbers, we noticed that

(1+0*I) / (+-Inf +- Inf*i) = (1+0*I) / (+-Inf + 0*i) = (1+0*I) / (0 +- Inf*i) = (NaN + NaN*I)

when using the type std::complexmpfr::mpreal. We think this is not desirable since, in these cases, some kind of "0" is warranted. The division also results in NaN if we replace (1+0*I) by (1+1*I) or by (0+0*I). For instance, using any of the built-in C++ datatypes std::complex, std::complex and std::complex, we obtain 0 on each of these operations. I am not saying the output 0 is preferable as opposed to NaN on all these cases, but I would expect at least (1+0*I)/(Inf+0*I) = 1/Inf = 0. [edited]

Is this the expected behavior when using std::complexmpfr::mpreal?
Is your intent to support std::complexmpfr::mpreal?

Other operations with undesired output are:

std::abs( (+-Inf +- Inf*i) ) = std::abs( (+-Inf +- Inf*i) ) = std::abs( (+-Inf +- Inf*i) ) = NaN,

where we expected the Inf propagation, and

std::abs( (0 + NaN*i) ) = 0,

where we expected the NaN propagation.

Thanks!

[advanpix]

Last time I checked implementation of std::complex<T> was very much dependent on 'double'-precision constants. Hence simple instantiation std::complex<mpfr::mpreal> would not work as expected (delivering low-accuracy results, etc.).
Basically we/you need to re-write overloaded class std::complex<mpfr::mpreal> from scratch.

As for the semantic of complex operations in degenerate cases (when go to infinity, division by zero, etc.) - this is (undefined) murky water. Every math package - MATLAB, Maple, Mathematica uses its own set of rules (contradicting each other). Add here different selection of branch cuts and you will see all hell breaks loose.

I would suggest to take a closer look on MPC library. Authors are doing a superb job on defining reasonable semantic for all these cases (but, again, it contradicts to other software packages:)).

I am closing the issue since it has nothing to do with the MPFR C++ wrapper on its own.

[weslleyspereira]

As for the semantic of complex operations in degenerate cases (when go to infinity, division by zero, etc.) - this is (undefined) murky water. Every math package - MATLAB, Maple, Mathematica uses its own set of rules (contradicting each other). Add here different selection of branch cuts and you will see all hell breaks loose.

I was aware of the contradictions between C and Fortran on Inf and NaN propagation. Thank you for highlighting those other examples.

I would suggest to take a closer look on MPC library. Authors are doing a superb job on defining reasonable semantic for all these cases (but, again, it contradicts to other software packages:)).

I will definitely take a look. Thanks!