Bug 221319

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.