Re: [PATCH] fortran/69121 -- Make IEEE_SCALB generic

[Steve Kargl <sgk@troutmask.apl.washington.edu>]

On Fri, Dec 21, 2018 at 11:07:08AM +0200, Janne Blomqvist wrote:
> On Fri, Dec 21, 2018 at 8:22 AM Steve Kargl <
> sgk@troutmask.apl.washington.edu> wrote:
> 
> > On Thu, Dec 20, 2018 at 01:47:39PM -0800, Steve Kargl wrote:
> > > The attached patch has been tested on x86_64-*-freebsd.
> > >
> > > OK to commit?
> > >
> > > 2018-12-20  Steven G. Kargl  <kargl@gcc.gnu.org>
> > >
> > >       PR fortran/69121
> > >       * libgfortran/ieee/ieee_arithmetic.F90: Provide missing functions
> > >       in interface for IEEE_SCALB.
> > >
> > > 2018-12-20  Steven G. Kargl  <kargl@gcc.gnu.org>
> > >
> > >       PR fortran/69121
> > >       * gfortran.dg/ieee/ieee_9.f90: New test.
> >
> > Now, tested on i586-*-freebsd.
> >
> 
> Hi, looks ok for trunk.
> 
> A few questions popped into my mind while looking into this:
> 
> 1) Why are none of the _gfortran_ieee_scalb_X_Y functions mentioned in
> gfortran.map? I guess they should all be there?
> 
> 2) Currently all the intrinsics map to the scalbn{,f,l} builtins. However,
> when the integer argument is of kind int64 or int128 we should instead use
> scalbln{,f,l}. This also applies to other intrinsics that use scalbn under
> the hood.
> 
> To clarify, fixing these is not a prerequisite for accepting the patch (I
> already accepted it), but more like topics for further work.
> 

I forgot to address your had a 2) item above.  ieee_scalb appears
to do the right thing.  FX addressed that with his implementation.
The 2nd argument is always cast to integer after reducing the range
to that of integer(4).

The binary floating point representation for a REAL(16) finite number
is x=f*2**e with f in [0.5,1) and e in [-16059,16384].  scalb(x,n) is
x*2**n, which becomes f*2**e*2**n = f*2**(e+n).  If x is the smallest
positive subnormal number, then n can be at most 32443 to still return 
a finite REAL(16) number.  Any larger value overflows to infinity.
If x is the largest positive finite number, then n can be -32443 to
return the small positive subnormal number.  Any more negative value
of n underflows to zero.  (Note, I could be off-by-one, but that is
just a detail.)

Consider

function foo(x,i)
   use ieee_arithmetic
   real(16) foo, c
   integer(8) i
   print *, ieee_scalb(c, i)
end function foo

-fdump-tree-original gives 
 
D.3853 = *i;
__result_foo = scalbnq (c,
    (integer(kind=4)) MAX_EXPR <MIN_EXPR <D.3853, 2147483647>, -2147483647>);

The range [-32443,32443] is a subset of [-huge(),huge(0)].

-- 
Steve