Division by complex<double> in clang++ versus g++

Question

When I compile the following code with g++ (4.8.1 or 4.9.0) or clang++ (3.4) I get different outputs.

#include <iostream>
#include <complex>

int main() {
  std::complex<double> c = {1.e-162,0};
  std::cout << 1.0/c << std::endl;
  return 0;
}

g++:

(1e+162,0)

clang++:

(inf,-nan)

Is this a bug in clang?

Update:

Thank you for your answers! I reported the bug: http://llvm.org/bugs/show_bug.cgi?id=19820

Answer 1

The standard says in [complex.numbers] (26.4/3):

If the result of a function is not mathematically defined or not in the range of representable values for its type, the behavior is undefined.

There are no specifics on how division should be implemented for complex numbers. Only in [complex.member.ops] it says:

complex<T>& operator/=(const complex<T>& rhs);

Effects: Divides the complex value rhs into the complex value *this and stores the quotient in *this. Returns: *this.

and in [complex.ops]:

template<class T> complex<T> operator/(const T& lhs, const complex<T>& rhs);

Returns: complex<T>(lhs) /= rhs.

---

As the inverse of 1.e-162 is 1.e+162 and this number is in the range of representable values for a double, the behavior is well defined.

Thus gcc gets it right and clang has a bug.

Answer 2

Ok so I have this wild guess, that in clang the complex number division is implemented like described on wiki: http://en.wikipedia.org/wiki/Complex_number#Multiplication_and_division.

One can see that the denominator is in the form c^2 + d^2. So 1.e-162 squared actually falls out of the IEE754 representable range for double which is std::numeric_limits<double>::min() - 2.22507e-308, and we have an underflow.

gcc somehow works this out, but if clang does simple square, as per @40two's standard quotation it enters into UB, and treats it as 0 after performing 1.e-162^2 + 0.0^2.

I tested clang and gcc for a number that should not result with underflow when squared.

#include <iostream>
#include <complex>

int main() {
  std::complex<double> c = {1.e-104,0};
  std::cout << 1.0/c << std::endl;
  return 0;
}

Results are fine:

luk32:~/projects/tests$ g++ --std=c++11 ./complex_div.cpp 
luk32:~/projects/tests$ ./a.out 
(1e+104,0)
luk32:~/projects/tests$ clang++ --std=c++11 ./complex_div.cpp 
luk32:~/projects/tests$ ./a.out 
(1e+104,0)

Not sure if this is a bug. But I think this is what is going on.

Addendum:

(inf,-nan) is also consistent if one evaluates those expressions by hand

We get:

real = (ac+bd) / (o)  - real part 
imag = (bc-ad) / (o)  - imaginary part

{a,b} = {1.0, 0.0}
{c,d} = {1.e-104, 0.0}

o is (c^2 + d^2) and we assume it is 0.0.

real / o = (1.e-104 * 1.0 + 0.0 * 0.0) / o = 1/o = inf
imag / o = (0.0 * 1.e-104 - 1.0 * 0.0) / o = -0.0 / o = -nan

I am just not absolutetly sure about the sign of -0.0 and -nan, I don't know IEE754 enough to evaluate (0.0 * 1.e-104 - 1.0 * 0.0). But everything seems consistent.