http://gcc.gnu.org/bugzilla/show_bug.cgi?id=56542
Bug #: 56542
Summary: complex number division underflow on darwin11 without
-lm
Classification: Unclassified
Product: gcc
Version: 4.2.1
Status: UNCONFIRMED
Severity: normal
Priority: P3
Component: c++
AssignedTo: unassig...@gcc.gnu.org
ReportedBy: pickledeyeb...@gmail.com
Created attachment 29589
--> http://gcc.gnu.org/bugzilla/attachment.cgi?id=29589
source code, *i file, screen output, compiler flags, etc.
This is similar to (closed) bug 42333 that was for darwin10 + -lm. This bug
report is for darwin11 without the math lib linked in and no apparent complex
number related compiler flags.
x/y where x, y < (1e-22,i*a) and x-y>eps yields (NaN,NaN) instead of 1.0.
numeric_limit<float>::eps() = 1.17549e-38
The essential operation of "complex::operator/" is like this:
(x + i*y) / (a + i*b)
= (x * i*y) * (a - i*b) / [(a + i*b) * (a - i*b)]
= [ (x*a + y*b) + (y*a - x*b) *i ] / (a*a + b*b)
= (x*a + y*b) / ((a*a * b*b) + i * (y*a - x*b) / (a*a + b*b)
It appears that the code is trying to rescale the problem, but the rescale is
failing due to underflow in denominator (and possibly numerator).
The following code demonstrates this bug:
# include <iostream>
# include <complex>
# include <cstdio>
# include <numeric>
using namespace std;
int main(void)
{
cout << "========================================================" << endl;
cout << "========================================================" << endl;
cout << "=== ===" << endl;
cout << "=== Demonstration of division bug in complex numbers ===" << endl;
cout << "=== ===" << endl;
cout << "========================================================" << endl;
cout << "========================================================" << endl;
cout << "limit<float>::min() = " << numeric_limits<float>::esp() << endl;
float base1 = 1.1754901;
float base2 = 1.1754900;
int largest_good = -22; // largest good exponent
int smallest_bad = -23; // smallest bad exponent
for(int exp=largest_good; exp>=smallest_bad; --exp) {
complex<float> x(base1*pow(10.,exp),0) ;
complex<float> y(base2*pow(10.,exp),0) ;
printf("Exponent: %d\n", exp);
printf("x = (%1.7fe%d, i*0)\n", base1, exp);
printf("y = (%1.7fe%d, i*0)\n", base2, exp);
printf("x / y = %1.16e / %1.16e\n", x.real(), y.real());
cout << "= " << x / y << endl;
cout << endl;
}
return 0;
}
It prints out:
========================================================
========================================================
=== ===
=== Demonstration of division bug in complex numbers ===
=== ===
========================================================
========================================================
limit<float>::eps() = 1.17549e-38
Exponent: -22
x = (1.1754901e-22, i*0)
y = (1.1754900e-22, i*0)
x / y = 1.1754900914587337e-22 / 1.1754899652409888e-22
= (1,0)
Exponent: -23
x = (1.1754901e-23, i*0)
y = (1.1754900e-23, i*0)
x / y = 1.1754901545676061e-23 / 1.1754899967954250e-23
= (nan,nan)
