Red Hat Bugzilla – Bug 221319

ieee arithmetics wrong

Last modified: 2007-11-30 17:11:52 EST

Description of problem: Devision of complex number by zero results in (nan, nan) instead of (inf, inf) Version-Release number of selected component (if applicable): gcc version 4.1.1 20061011 (Red Hat 4.1.1-30) How reproducible: Always Steps to Reproduce: 1. Compile the code below 2. run 3. Actual results: (nan, nan) Expected results: (inf, inf) Additional info: This is a test code: #include <complex> #include <iostream> int main (void) { std::cout << std::complex<double> (1.0, 1.0) / 0.0 << std::endl; return 0; }

I skimmed the C++ standard and didn't find where it would require particular behavior for this. ISO C99 specifies this (provides exact algorithm for when one or both operands of * resp. / are not complex and additionally requires that division by 0 + 0 * I results in infinity (if the first operand is non-zero finite or infinity), though, infinity in ISO C99 is anything where at least one of the real and imaginary parts are inf. There is a bug in libgcc's __div[sdtx]c3 (PR30360), but that's not what matters for this testcase. The testcase gives (inf, inf) when compiled with optimizations (because there the optimizer detects this division has a real second operand rather than complex and can use the cheaper (a, b) / (c, 0) = (a / c, b / c) algorithm), but at -O0 find_lattice_value is not called (for simple cases like this it wouldn't be hard to detect it even at -O0 though) and full complex / complex division is used. As C++ is not using flag_complex_method = 2 (which gcc -std=c99 or -std=gnu99 uses), even when libgcc is fixed this will still give NaN. So, before doing anything further it would be good to understand what exactly ISO C++ requires here.

E.g. if ISO C++ required that non-zero finite or infinite complex<float> when divided by float 0.0 (i.e. when the operator/= (float) is used; similarly for double and long double) gives infinity, but no requirements about the other operator (with complex<float> argument), then we can either change tree-complex.c to do some extra work to find ONLY_REAL/ONLY_IMAG/etc. at -O0, or we can help it on the libstdc++-v3 headers side, like this: --- libstdc++-v3/include/std/complex 2006-12-08 15:58:11.000000000 +0100 +++ libstdc++-v3/include/std/complex 2007-01-04 11:38:14.000000000 +0100 @@ -1097,7 +1097,12 @@ _GLIBCXX_BEGIN_NAMESPACE(std) inline complex<float>& complex<float>::operator/=(float __f) { - _M_value /= __f; + // _M_value /= __f; + // Help GCC here a little bit, when not optimizing + // it wouldn't find that the second argument is not complex. + _ComplexT __t = _M_value; + __real__ _M_value = __real__ __t / __f; + __imag__ _M_value = __imag__ __t / __f; return *this; } @@ -1250,7 +1255,12 @@ _GLIBCXX_BEGIN_NAMESPACE(std) inline complex<double>& complex<double>::operator/=(double __d) { - _M_value /= __d; + // _M_value /= __d; + // Help GCC here a little bit, when not optimizing + // it wouldn't find that the second argument is not complex. + _ComplexT __t = _M_value; + __real__ _M_value = __real__ __t / __d; + __imag__ _M_value = __imag__ __t / __d; return *this; } @@ -1403,7 +1413,12 @@ _GLIBCXX_BEGIN_NAMESPACE(std) inline complex<long double>& complex<long double>::operator/=(long double __r) { - _M_value /= __r; + // _M_value /= __r; + // Help GCC here a little bit, when not optimizing + // it wouldn't find that the second argument is not complex. + _ComplexT __t = _M_value; + __real__ _M_value = __real__ __t / __r; + __imag__ _M_value = __imag__ __t / __r; return *this; } Or, if ISO C++ requires that even complex<double> (1.0, 1.0) / complex<double> (0.0, 0.0) gives infinity, we'd need to use flag_complex_method = 2 even in ISO C++.

The issue here is not C++ standard compliance, but rather IEEE 754 standard complience. Does C++ (or C99) standard codify what the result of operations like x = 3.0/2.0 should be? Sincerely, Dmitri.

Per the C++ Standard, division by zero is undefined [expr]/5, but it says there may be library functions that can control the results. I suppose this is meant to defer to other standards. [lib.complex.member.ops]/7 and /15 define division of complex by scalar and by complex, respectively, without any specifics as to the arithmetic approach, which further supports this theory.

I do not have access to C++ language standard, but I guess the language standard cannot require IEEE floating math compliance. But I would expect that language _implementation_ on IEEE 754 compliant hardware to be complient as well. Again, what does the C++ standard say about accuracy of floating point arithmetic? Sincerely, Dmitri.

Moving further discussions upstream.