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

[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))'.

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

[rguenth at gcc dot gnu.org]

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

--- Comment #1 from Richard Guenther <rguenth at gcc dot gnu.org> 2011-04-23 09:20:46 UTC ---
Integer exponent powers are computed using

TYPE
NAME (TYPE x, int m)
{
  unsigned int n = m < 0 ? -m : m;
  TYPE y = n % 2 ? x : 1;
  while (n >>= 1)
    {
      x = x * x;
      if (n % 2)
        y = y * x;
    }
  return m < 0 ? 1/y : y;
}

or optimal power series expansion for constant exponents starting with
GCC 4.5.

ISTR the standard allows this.  You should use a double exponent to
use the real power function, pow (2., -1024.)

[Bug libstdc++/48738] 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

--- Comment #2 from Paulo A. P. Pires <paulo1205 at yahoo dot com> 2011-04-23 13:55:48 UTC ---
(In reply to comment #1)
> Integer exponent powers are computed using
> 
> TYPE
> NAME (TYPE x, int m)
> {
>   unsigned int n = m < 0 ? -m : m;
>   TYPE y = n % 2 ? x : 1;
>   while (n >>= 1)
>     {
>       x = x * x;
>       if (n % 2)
>         y = y * x;
>     }
>   return m < 0 ? 1/y : y;
> }
> 
> or optimal power series expansion for constant exponents starting with
> GCC 4.5.
> 
> ISTR the standard allows this.  You should use a double exponent to
> use the real power function, pow (2., -1024.)

I have no objection to using such optimizations for integer exponents, and I am
aware that it .  But perhaps it could be improved further, by widening the
range over which the result for an integer exponent will be consistent with the
corresponding double (or whatever) exponent.

If the code above was changed to something like what goes below, would it be
worse?

TYPE
NAME (TYPE x, int m)
{
  unsigned int n;
  if(m < 0){
    n = -m;
    x = TYPE(1) / x;
  }
  TYPE y = n % 2 ? x : 1;
  while (n >>= 1)
    {
      x = x * x;
      if (n % 2)
        y = y * x;
    }
  return y;
}

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

[rguenth at gcc dot gnu.org]

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

--- Comment #3 from Richard Guenther <rguenth at gcc dot gnu.org> 2011-04-23 16:24:08 UTC ---
(In reply to comment #2)
> (In reply to comment #1)
> > Integer exponent powers are computed using
> > 
> > TYPE
> > NAME (TYPE x, int m)
> > {
> >   unsigned int n = m < 0 ? -m : m;
> >   TYPE y = n % 2 ? x : 1;
> >   while (n >>= 1)
> >     {
> >       x = x * x;
> >       if (n % 2)
> >         y = y * x;
> >     }
> >   return m < 0 ? 1/y : y;
> > }
> > 
> > or optimal power series expansion for constant exponents starting with
> > GCC 4.5.
> > 
> > ISTR the standard allows this.  You should use a double exponent to
> > use the real power function, pow (2., -1024.)
> 
> I have no objection to using such optimizations for integer exponents, and I am
> aware that it .  But perhaps it could be improved further, by widening the
> range over which the result for an integer exponent will be consistent with the
> corresponding double (or whatever) exponent.
> 
> If the code above was changed to something like what goes below, would it be
> worse?
> 
> TYPE
> NAME (TYPE x, int m)
> {
>   unsigned int n;
>   if(m < 0){
>     n = -m;
>     x = TYPE(1) / x;
>   }
>   TYPE y = n % 2 ? x : 1;
>   while (n >>= 1)
>     {
>       x = x * x;
>       if (n % 2)
>         y = y * x;
>     }
>   return y;
> }

I'm sure you can create a similar issue for this variant.  Fact is that
there are multiple roundings involved in computing integer powers.

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

[paolo.carlini at oracle dot com]

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

--- Comment #4 from Paolo Carlini <paolo.carlini at oracle dot com> 2011-04-23 20:26:47 UTC ---
I'm following with interest the numerical analysis side of the discussion, but
I have to add that, *as a libstc++-v3 proper* issue, as mentioned by Richard,
the issue is moot, because in currently active release branches either the
library uses a front-end builtin (__builtin_powi*) or just pow unconditionally
(in C++0x mode). In other terms, *minimal* open-coding in the library for these
functions (and we mean to stick to this general design choice, AFAIK)

[Bug tree-optimization/48738] pow() fails to produce (some) subnormalized numbers with integer exponents

[paolo.carlini at oracle dot com]

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

Paolo Carlini <paolo.carlini at oracle dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|paulo1205 at yahoo dot com  |paolo.carlini at oracle dot
                   |                            |com, wschmidt at gcc dot
                   |                            |gnu.org
          Component|libstdc++                   |tree-optimization

--- Comment #5 from Paolo Carlini <paolo.carlini at oracle dot com> 2011-05-25 09:54:15 UTC ---
I don't think this is a library issue. Maybe Bill is interested / has an
opinion..

[Bug tree-optimization/48738] pow() fails to produce (some) subnormalized numbers with integer exponents

[wschmidt at gcc dot gnu.org]

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

--- Comment #6 from William J. Schmidt <wschmidt at gcc dot gnu.org> 2011-05-25 12:15:46 UTC ---
(In reply to comment #5)
> I don't think this is a library issue. Maybe Bill is interested / has an
> opinion..

Just to agree that it is WAD.  The proposed solution just moves the values of x
and m for which the rounding will cause edge-case anomalies.  Either way, the
compounded rounding error means the two versions of pow will never produce
identical results.  When precision is important, use the real-valued version.