13

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

Michel Ferrero
  • 131
  • 5
  • You should check whether values are within the numeric limits, just to be on the safe side. – luk32 May 07 '14 at 13:39
  • I checked and they are, the doubles range from 2.22507e-308 to 1.79769e+308 – Michel Ferrero May 07 '14 at 13:46
  • I checked it too, std::defaultfloat << std::numeric_limits::min() gives 2.22507e-308 so this should not be an issue. If it went into a denormalized number it could behave like that. – luk32 May 07 '14 at 13:47
  • I smell an IEEE 754 issue... and underflow – Mgetz May 07 '14 at 13:47

4 Answers4

5

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.

Danvil
  • 22,240
  • 19
  • 65
  • 88
4

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.

luk32
  • 15,812
  • 38
  • 62
2

Update:

Quoting from the standard:

If during the evaluation of an expression, the result is not mathematically defined or not in the range of representable values for its type, the behavior is undefined. [ Note: most existing implementations of C++ ignore integer overflows. Treatment of division by zero, forming a remainder using a zero divisor, and all floating point exceptions vary among machines, and is usually adjustable by a library function].

Quoting from http://lists.freebsd.org/pipermail/freebsd-numerics/2014-March/000549.html:

It appears that clang developers have chosen the naive complex division algorithm.

...

I did a bit of grepping. Could it be that the division algorithm is contained in the file src/contrib/llvm/tools/clang/lib/CodeGen/CGExprComplex.cpp inside the function ComplexExprEmitter::EmitBinDiv ?

If you look at the code, it certainly looks like it is generating code to perform complex division, and it definitely looks like they are using the naive algorithm.

Assuming that indeed the clang uses naive complex division the expression 1.0 / c evaluates according to the naive implementation of complex division to the following expression enter image description here,

enter image description here

1.e-324 is out of the double range. This results according to the standard to undefined behaviour.

Also making a search in the LLVM/Clang bug list, it appears that there are quite some issues concerning complex division.

As such your case is a bug and you should report it.

For anyone who is interested on how robust complex division is implemented take a look at

  1. http://ideone.com/bqFk8j and
  2. A Robust Complex Division in Scilab.
101010
  • 41,839
  • 11
  • 94
  • 168
  • 1
    Why is this not mathematically defined or out of range? – Mat May 07 '14 at 13:40
  • I think the `1e-162` is smaller then smallest possible representation in double. Thus, the compiler is free to do whatever on the division. gcc is nice, and clang treats it as `0`. But according to the quote it is allowed to. The reason why gcc works with it at all is because it is a [denormalized number](http://en.wikipedia.org/wiki/Denormal_number) – luk32 May 07 '14 at 13:41
  • @luk32: "The 11 bit width of the exponent allows the representation of numbers with a decimal exponent between 10^−308 and 10^308, with full 15-17 decimal digits precision." – Danvil May 07 '14 at 13:42
  • @Danvil OK, `std::defaultfloat << std::numeric_limits::min()` gives `2.22507e-308` so this should not be an issue. – luk32 May 07 '14 at 13:46
  • 2
    Down-voters please explain the reason of your down-vote. – 101010 May 07 '14 at 14:20
  • `a / b` is defined as the number `x` for which we have `a = b * x`. Thus for `a=1` and `b=1e-162`, the result is `1e+162` which is in the range of representable values. – Danvil May 07 '14 at 14:26
  • Well I think, that pure quotation of standard did not explain anything. The other answers suggest there is nothing how the division should be implemented. It only says that if the **result** result is not representable it is UB. And the result of this division should be representable. So it looks like they implemented it "naively" with square operation, so it fails sometimes. The same goes for `1e+162`. `clang` gives `(0,0)`, while gcc `(1e-162,0)`. I think it's somewhat justified to call it implementation bug. – luk32 May 07 '14 at 14:27
  • @40two regarding http://ideone.com/ClGxhp, assuming that you are the author, could you please correct the test to include both real and imag, at the moment it compares real(b)>=real(b). – Lutz Lehmann May 08 '14 at 16:43
  • The standard text about UB refers to stuff in the user code. It doesn't apply to the implementation itself. (I agree that it is a compiler bug, but because "the implementation is wrong", not "the implementation causes UB"). – M.M May 21 '14 at 23:19
1

Yes. You should report a bug.

The complex number 1e-162+0i is identical to the real number 1e-162. This falls well within the range of double. The reciprocal of this value is and should be 1e+162. Any other value represents a bug in the arithmetic libraries.

The value returned by clang is wrong. It's a bug.

david.pfx
  • 10,520
  • 3
  • 30
  • 63