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

["rguenth at gcc dot gnu.org" <gcc-bugzilla at gcc dot gnu dot 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.