fortran/69121 -- Make IEEE_SCALB generic

The attached patch has been tested on x86_64-*-freebsd.

It provides the missing pieces for IEEE_SCALB to allow
it to be generic function.  Briefly, when FX contributed
IEEE_SCALB, he implemented a casting of the integer
argument to the default integer kind.  For integer(1)
and integer(2), this is a simple cast.  For integer(8)
and integer(16), he first checks to see if the value 
is in the range of default integer kind.  So, one
essentially has

ieee_scalb(x, 2_1) --> ieee_scalb(x, int(2_1, 4)) 
ieee_scalb(x, 2_8) --> ieee_scalb(x, int(max(min(2_8,huge(0)),-huge(0)), 4)

Unfortnately, the interface exposed by the IEEE_ARITHMETIC
module did not allow for non-default integer kind arguments.
The patch corrects that omission.

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.

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.

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.
>
> --
> Steve
>

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.

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.

Just, sent shorter note in private email to Thomas.

No, I'm adding the missing functions to the INTERFACE.

This will not compile: 

   program foo
   use ieee_arithmetic
   real x
   integer(8) i
   x = 2
   i = 2_8
   print *, ieee_scalb(x,i)
   end program

because the module has a generic interface that does not
include the integer(8) argument.

FX seems to have wanted to avoid the explosion of functions in
the library.  In trans-intrinsic.c (conv_intrinsic_ieee_scalb),
he does a conversion of the integer(8) to an integer(4).
Unfortunately, checking the interface for ieee_scalb occurs
before code generation.  Compiling the above after my patch,
the -ftree-dump-original contains

D.3769 = __builtin_scalbnf (x, (integer(kind=4)) MAX_EXPR <MIN_EXPR <D.3768, 2147483647>, -2147483647>);

Hi Steve,

> No, I'm adding the missing functions to the INTERFACE.

Ah, I see. What I missed was that the function is actually translated
to something else.

So, OK for trunk, and thanks for the patch!

Regards

	Thomas

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)].

On Fri, Dec 21, 2018 at 06:31:27PM +0100, Thomas Koenig wrote:
> Hi Steve,
> 
> > No, I'm adding the missing functions to the INTERFACE.
> 
> Ah, I see. What I missed was that the function is actually translated
> to something else.
> 
> So, OK for trunk, and thanks for the patch!
> 

Janne, Thomas, Thanks for the quick review.

I first reported the bug 2 years ago (2016-01-02).  I left
it as low-hanging fruit with the hope that a junior gfortran
hacker would cut his/her teeth on a patch.  Junior hasn't
come along, and with nearly 1000 open PRs, I decided to
fix it.  I took a 5-6 month break from working on gfortran.
When I looked a month ago, there were 1002 open PRs.  With my
recent burst of fixing/closing 16 PRs in December, we're 
still at 991.  Too many PRs, too few hands. :-(

On Fri, Dec 21, 2018 at 7:59 PM Steve Kargl <
sgk@troutmask.apl.washington.edu> wrote:

> 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)].
>

True. I guess the advantage of scalbln* would be to avoid the MAX/MIN_EXPR
and casting for kind int64.

On Fri, Dec 21, 2018 at 08:59:04PM +0200, Janne Blomqvist wrote:
> On Fri, Dec 21, 2018 at 7:59 PM Steve Kargl <
> >
> > 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)].
> >
> 
> True. I guess the advantage of scalbln* would be to avoid the MAX/MIN_EXPR
> and casting for kind int64.
> 

fdlibm-based libm has 

#define NMAX    65536
#define NMIN    -65536
double
scalbln(double x, long n)
{

        return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}

A search for glibc's libm locates https://tinyurl.com/ybcy8w4t
which is a bit-twiddling routine.  Not sure it's worth the
effort.  Joseph Myers might have an opinion.

On Fri, 21 Dec 2018, Steve Kargl wrote:

> scalbln(double x, long n)
> {
> 
>         return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
> }
> 
> A search for glibc's libm locates https://tinyurl.com/ybcy8w4t
> which is a bit-twiddling routine.  Not sure it's worth the
> effort.  Joseph Myers might have an opinion.

Such comparisons are needed in the scalbn / scalbln implementations anyway 
to deal with large exponents.  I suppose where there's a suitable scalbln 
implementation, and you don't know if the arguments are within the range 
of int, calling scalbln at least saves code size in the caller and avoids 
duplicating those range checks.

On Fri, Dec 21, 2018 at 07:39:45PM +0000, Joseph Myers wrote:
> On Fri, 21 Dec 2018, Steve Kargl wrote:
> 
> > scalbln(double x, long n)
> > {
> > 
> >         return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
> > }
> > 
> > A search for glibc's libm locates https://tinyurl.com/ybcy8w4t
> > which is a bit-twiddling routine.  Not sure it's worth the
> > effort.  Joseph Myers might have an opinion.
> 
> Such comparisons are needed in the scalbn / scalbln implementations anyway 
> to deal with large exponents.  I suppose where there's a suitable scalbln 
> implementation, and you don't know if the arguments are within the range 
> of int, calling scalbln at least saves code size in the caller and avoids 
> duplicating those range checks.
> 

I was thinking along the lines of -ffast-math and whether 
__builtin_scalbn and __builtin_scalbln are then inlined.  
The comparisons may inhibit inlining __builtin_scalbn;
while, if gfortran used __builtin_scalbln, inlining would
occur.

As it is, for

   function foo(x,i)
     use ieee_arithmetic
     real(8) foo, c
     integer(8) i
     foo = ieee_scalb(c, i)
   end function foo

the options -ffast-math -O3 -fdump-tree-optimized give

  <bb 2> [local count: 1073741824]:
  _gfortran_ieee_procedure_entry (&fpstate.0);
  _8 = *i_7(D);
  _1 = MIN_EXPR <_8, 2147483647>;
  _2 = MAX_EXPR <_1, -2147483647>;
  _3 = (integer(kind=4)) _2;
  _4 = __builtin_scalbn (c_9(D), _3);
  _gfortran_ieee_procedure_exit (&fpstate.0);
  fpstate.0 ={v} {CLOBBER};
  return _4;

It seems this could be 

  <bb 2> [local count: 1073741824]:
  _gfortran_ieee_procedure_entry (&fpstate.0);
  _3 = (integer(kind=4)) *i_7(D);
  _4 = __builtin_scalbn (c_9(D), _3);
  _gfortran_ieee_procedure_exit (&fpstate.0);
  fpstate.0 ={v} {CLOBBER};

Hi Steve

On 21/12/18 8:01 PM, Steve Kargl wrote:
> On Fri, Dec 21, 2018 at 07:39:45PM +0000, Joseph Myers wrote:
>> On Fri, 21 Dec 2018, Steve Kargl wrote:
>>
>>> scalbln(double x, long n)
>>> {
>>>
>>>          return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
>>> }
>>>
>>> A search for glibc's libm locateshttps://tinyurl.com/ybcy8w4t
>>> which is a bit-twiddling routine.  Not sure it's worth the
>>> effort.  Joseph Myers might have an opinion.
>> Such comparisons are needed in the scalbn / scalbln implementations anyway
>> to deal with large exponents.  I suppose where there's a suitable scalbln
>> implementation, and you don't know if the arguments are within the range
>> of int, calling scalbln at least saves code size in the caller and avoids
>> duplicating those range checks.
>>
> I was thinking along the lines of -ffast-math and whether
> __builtin_scalbn and __builtin_scalbln are then inlined.
> The comparisons may inhibit inlining __builtin_scalbn;
> while, if gfortran used __builtin_scalbln, inlining would
> occur.
>
> As it is, for
>
>     function foo(x,i)
>       use ieee_arithmetic
>       real(8) foo, c
>       integer(8) i
>       foo = ieee_scalb(c, i)
>     end function foo
>
> the options -ffast-math -O3 -fdump-tree-optimized give
>
>    <bb 2> [local count: 1073741824]:
>    _gfortran_ieee_procedure_entry (&fpstate.0);
>    _8 = *i_7(D);
>    _1 = MIN_EXPR <_8, 2147483647>;
>    _2 = MAX_EXPR <_1, -2147483647>;
>    _3 = (integer(kind=4)) _2;
>    _4 = __builtin_scalbn (c_9(D), _3);
>    _gfortran_ieee_procedure_exit (&fpstate.0);
>    fpstate.0 ={v} {CLOBBER};
>    return _4;
>
> It seems this could be
>
>    <bb 2> [local count: 1073741824]:
>    _gfortran_ieee_procedure_entry (&fpstate.0);
>    _3 = (integer(kind=4)) *i_7(D);
>    _4 = __builtin_scalbn (c_9(D), _3
>    _gfortran_ieee_procedure_exit (&fpstate.0);
>    fpstate.0 ={v} {CLOBBER};
>
I am observing your new test pr88328.f90 failing on 
arm-none-linux-gnueabihf:
Excess errors:
/build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:20:36: Error: 
Invalid kind for REAL at (1)
/build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:35:36: Error: 
Invalid kind for REAL at (1)
/build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:50:36: Error: 
Invalid kind for REAL at (1)
/build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:65:36: Error: 
Invalid kind for REAL at (1)

On Mon, Dec 24, 2018 at 12:26:46PM +0000, Sudakshina Das wrote:
> >
> I am observing your new test pr88328.f90 failing on 
> arm-none-linux-gnueabihf:
> Excess errors:
> /build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:20:36: Error: 
> Invalid kind for REAL at (1)
> /build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:35:36: Error: 
> Invalid kind for REAL at (1)
> /build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:50:36: Error: 
> Invalid kind for REAL at (1)
> /build/src/gcc/gcc/testsuite/gfortran.dg/ieee/ieee_9.f90:65:36: Error: 
> Invalid kind for REAL at (1)

pr88328.f90 or ieee_9.f90?

In either case, you can XFAIL the testcase for arm-none-linux-gnueabihf.