[Bug c++/56542] New: complex number division underflow on darwin11 without -lm

["pickledeyeball at gmail dot com" <gcc-bugzilla at gcc dot gnu dot org>]

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: unassigned@gcc.gnu.org
        ReportedBy: pickledeyeball@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)