Hello I have found a bug in the implementation of the complex library of g++ and the complex.h library of the gcc (c99 support).
The simplest program that demonstrates the bug is: -------------------------------- #include<iostream> #include<complex> using namespace std; int main() { complex<double> a = 1.0/0.0; // + infinity complex<double> b = 1.0; complex<double> c = b/a; // should be 0 but is (NaN,NaN) cout << a << endl; cout << b << endl; cout << c << endl; } -------------------------------- Since all values are real one should expect the result 0 (the IEEE floating-point value for 1/infinity). More serious is that if a contains a large but representable value for which |a|^2 is not representable the implementation will yields bad results. For example let a be the largest representable number and b be a/2 then b/a should be 0.5 and not NaN. The bug comes from the naive implementation a + b I ac + bd bc - ad ------- = ------- + ------- I c + d I c^2+d^2 c^2+d^2 We get unnecessary overflow errors if c^2+d^2 is too large. For a better result one should use the following formula due to Smith (see "What Every Computer Scientist Should Know About Floating-Point Arithmetic" http://cch.loria.fr/documentation/IEEE754/) | a + b(d/c) c - a(d/c) | ---------- + ---------- I if |d| < |c| a + b I | c + d(d/c) a + d(d/c) ------- = | c + d I | b + a(c/d) -a+ b(c/d) | ---------- + ---------- I if |d| >= |c| | d + c(c/d) d + c(c/d) Even better we should test for |d| = |c| and definie a + b I a + b s c - a s ------- = ------- + ------- if |d| = |c| c + d I c + d s d + c s where s = 1 if c and d have the same sign bit and s = -1 otherwise. (Especially if c = +0.0 and d = -0.0 we must set s = -1.) This improves the results if |c|=|d|=0 or |c|=|d|=infinity. So my suggestion is to add template specialization for the operator/= for the types complex<float> complex<double> and complex<long double> that use Smith formula. (For integer types I think we should stay with the old implementation.) I would volunteer to fix the c++ library. But I have less experience in C99 so it would be better if another can fix the c library. I guess that the same bug will be present in the libraries of the other languages (Fortan, etc.), but I have not tested it. I am sorry but I have difficulties with the bug database, so I send you the mail directly. With best regards Andreas Klein Institute for Mathematics and Computer Science [EMAIL PROTECTED] University of Kassel Germany