complex numbers

Andreas Klein klein@mathematik.uni-kassel.de
Mon Dec 6 09:11:00 GMT 2004


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
    klein@mathematik.uni-kassel.de
    University of Kassel
    Germany



More information about the Gcc-bugs mailing list