Bug 221319 - ieee arithmetics wrong
ieee arithmetics wrong
Product: Fedora
Classification: Fedora
Component: gcc (Show other bugs)
All Linux
medium Severity medium
: ---
: ---
Assigned To: Jakub Jelinek
Depends On:
  Show dependency treegraph
Reported: 2007-01-03 14:47 EST by Dmitri A. Sergatskov
Modified: 2007-11-30 17:11 EST (History)
3 users (show)

See Also:
Fixed In Version:
Doc Type: Bug Fix
Doc Text:
Story Points: ---
Clone Of:
Last Closed: 2007-01-16 05:23:43 EST
Type: ---
Regression: ---
Mount Type: ---
Documentation: ---
Verified Versions:
Category: ---
oVirt Team: ---
RHEL 7.3 requirements from Atomic Host:
Cloudforms Team: ---

Attachments (Terms of Use)

External Trackers
Tracker ID Priority Status Summary Last Updated
GNU Compiler Collection 30482 None None None Never

  None (edit)
Description Dmitri A. Sergatskov 2007-01-03 14:47:57 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:

Steps to Reproduce:
1. Compile the code below
2. run
Actual results:
(nan, nan)

Expected results:
(inf, inf)

Additional info:

This is a test code:
 #include <complex>
 #include <iostream>

 main (void)
   std::cout << std::complex<double> (1.0, 1.0) / 0.0 << std::endl;
   return 0;
Comment 1 Jakub Jelinek 2007-01-04 05:30:52 EST
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.
Comment 2 Jakub Jelinek 2007-01-04 05:44:58 EST
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++.
Comment 3 Dmitri A. Sergatskov 2007-01-04 10:52:18 EST
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?

Comment 4 Alexandre Oliva 2007-01-04 14:30:24 EST
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.
Comment 5 Dmitri A. Sergatskov 2007-01-04 18:13:23 EST
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?


Comment 6 Jakub Jelinek 2007-01-16 05:23:43 EST
Moving further discussions upstream.

Note You need to log in before you can comment on or make changes to this bug.