[Bug libstdc++/48738] New: pow() fails to produce (some) subnormalized numbers with integer exponents

[paulo1205 at yahoo dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48738

           Summary: pow() fails to produce (some) subnormalized numbers
                    with integer exponents
           Product: gcc
           Version: 4.4.3
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: paulo1205@yahoo.com


The standard <cmath>'s function pow() fails to produce some sub-normalized
numbers if the version that takes an integer exponent is used (using G++ 4.4.3,
x86_64-linux-gnu).

For example, the following piece of code will cause 'p' to get a zero value,
instead of the expected sub-normalized value 5.56268464626e-309.

    double p=pow(2.0, -1024);

I dug the include files to try to discover the reason for this strange result. 
It seemed to me that the function that gets called when the exponent is integer
first checks if the exponent is negative, and calls yet another overloaded
version that takes an unsigned exponent, passing to it the absolute value of
the original exponent, and dividing 1.0 by the received result.

So, apparently my original call 'pow(2.0, -1024)' became translated to
'1.0/pow(2.0, 1024u)'.  However, as 2**1024 is too large for a double, it
becomes the positive infinity, and the result of dividing 1.0 by this infinity
is zero.

As a sanity check, I tried to make the exponent a double (i.e. I did

    double p=pow(2.0, -1024.0);

instead).  Not surprisingly, now I got 5.56268464626e-309, which was the result
that I expected from the beginning with an integer exponent.



POSSIBLE FIX:

Maybe the fix is simple (yet I didn't work it to a very long extent): instead
of making 'pow(some_double, some_negative_int)' become '1.0/pow(some_double,
unsigned(-some_negative_int))', change it to 'pow(1.0/some_double,
unsigned(-some_negative_int))'.