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