Re: complex numbers

[Gabriel Dos Reis <gdr at integrable-solutions dot net>]

Andreas Klein <klein@mathematik.uni-kassel.de> writes:

| 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.
| 

As a matter of fact, the implementation of <complex> is criticized,
once in a while, because it does NOT use the grammar school rule you
present above.  However, for float, double, long double it specializes
to __complex__ T which is  what the compiler uses to implement complex
numbers in C99 and Fortran.  So, the problem is a compiler problem not
libstdc++ problem.

[...]

| 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.)

Did you look at the actual implementation?

-- Gaby