[J3] [BULK] SCALE(x,i) with complex x

Tue Nov 7 06:32:44 UTC 2023

On Tue, Nov 07, 2023 at 12:39:56PM +0900, Malcolm Cohen wrote:
> >2.**n * (1., 1.) will raise FE_INEXACT.
> 
> Not so. Every intermediate calculation there is exact.
> 

Whoops.  I should have been more careful.

> call test(3)

  call test(-128) // nagfor results?

> end
> subroutine test(n)
> use ieee_arithmetic
> complex c,d
> logical inex
> c = (1,1)
> d = 2.0**n*c
> call ieee_get_flag(ieee_inexact,inex)
> print *,inex,d
> end
> 
> This prints "f (8.000,8.000)"  (with variable numbers of zeroes) on gfortran, intel fortran, nag fortran, ...
> 
> >2.**n * (1., inf) = (nan,inf) will raise FE_INVALID.
> 
> Highly unlikely unless you use a pessimising compiler. Wow, gfortran is
> a pessimising compiler here?! I am very surprised.

"Pessimising" might be a bit strong.  gfortran implements real*complex
multiplication according the Fortran standard and correctly.  Type
conversion gives 2.**n*(x,y) -> (2.**n, 0.) * (x,y).  Short-circuiting
2.**n*(x,y) to (2.**n*x, 2.**n*y) by other compilers seems to suggest
those compilers are appealing to 10.1.5.2.4 where one is dismissing
the possibility that either x or y is exceptional.  Or, those 
compilers are using a distributive property of ordinary math.  Seems
to be QoI.

BTW, you get the NAG and Intel results with gfortran if one
uses the -ffast-math option.

> SCALE is slower than multiplication (unless of course, the compiler
> can see all the values and turn it into a multiply, in which case
> it's merely more opaque). Hardware multiplication is highly optimised.
> We are no longer in the days when cpus sometimes did not have hardware
> multiplication.

It's not that much slower.  But it can be more robust.
What does NAG and Intel do with 2**n*(x,y) if n > maxexponent(x)
and x and y are sufficiently less than 1.?

2.**(maxexponent(t0) + 1) * (0.125, 0.125)

I suppose forcing a user to do 2.**maxexponent(t0) * (0.125, 0.125) * 2
is acceptable.

> 
> Cheers,
> -- 
> ..............Malcolm Cohen, NAG Oxford/Tokyo.
> 
> -----Original Message-----
> From: J3 <j3-bounces at mailman.j3-fortran.org> On Behalf Of Steven G. Kargl via J3
> Sent: Tuesday, November 7, 2023 4:44 AM
> To: Long, Bill F <william.long at hpe.com>
> Cc: Steven G. Kargl <kargl at uw.edu>; General J3 interest list <j3 at mailman.j3-fortran.org>
> Subject: Re: [J3] [BULK] SCALE(x,i) with complex x
> 
> 2.**n * (1., 1.) will raise FE_INEXACT.
> 2.**n * (1., inf) = (nan,inf) will raise FE_INVALID.
> 
> SCALE((1.,1.), 2) is exact, .i.e, no exception raised.
> SCALE((1.,inf), 2) = (4., inf) will raise FE_OVERFLOW.
> 
> SCALE((0.125, 0.125), maxexponent(t0) + 1) = (8.507059173E+37,8.507059173E+37)
> 
> 2.**(maxexponent(t0) + 1) * (0.125, 0.125) does not compile.
> 
> My actual use case is computation of bessel functions
> of a complex argument via downward recursion.  Rescaling
> the recursion sequence allows one to extend the length
> of the sequence while avoiding overflow.
> 
> And, yes, I realize that this would be for F202Y.
> 
> -- 
> steve
> 
> On Mon, Nov 06, 2023 at 04:21:11PM +0000, Long, Bill F wrote:
> > In F90 there was a category of intrinsics "Floating point manipulation functions" that included EXPONENT, FRACTION, NEAREST, RRSPACING, SCALE, SET_EXPONENT, and SPAING.
> > For most of these only REAL arguments make sense. For SCALE, there is also a sensible result definition for a COMPLEX argument.  But given the context, I suspect that was not
> > in the scope of the conversation.  This could be added to the requirements list for F202Y.  (It is too late in the process to add this to F2023.)  But the question for F202Y would be "Why bother?".
> > 
> > Cheers,
> > Bill
> > 
> > ________________________________
> > From: J3 <j3-bounces at mailman.j3-fortran.org> on behalf of Clune, Thomas L. (GSFC-6101) via J3 <j3 at mailman.j3-fortran.org>
> > Sent: Monday, November 6, 2023 6:55 AM
> > To: kargl at uw.edu <kargl at uw.edu>; General J3 interest list <j3 at mailman.j3-fortran.org>
> > Cc: Clune, Thomas L. (GSFC-6101) <thomas.l.clune at nasa.gov>
> > Subject: Re: [J3] [EXTERNAL] [BULK] SCALE(x,i) with complex x
> > 
> > 
> > Hi Steve,
> > 
> > 
> > 
> > Extending to complex would not seem to pose any challenge.    As to why (well before my time), my guess would be that someone was attempting to emulate a specific bit of functionality in another layer (IEEE?) that did not address the complex case.
> > 
> > 
> > 
> > Of course back in the day, SCALE() could provide some significant performance advantage via bit-twiddling vs full multiply.     These days, though, unless an arg and result are known to be very small and cache resident, the memory access will dominate.   Given that I think a variant of your 3rd version is actually preferable in most implementations:
> > 
> > 
> > 
> > PROGRAM FOO
> > 
> >    COMPLEX A, B
> > 
> >    A = (1,1)
> > 
> >    B = (2**2) * A   ! parens added for clarity for reader
> > 
> >    PRINT *, B
> > 
> > END PROGRAM FOO
> > 
> > 
> > 
> > I find this implementation much more obvious than SCALE() as one needs to 1st remember which argument is which.   Of course, if I used SCALE() frequently, I might feel differently.
> > 
> > 
> > 
> > Not opposed to the extension – just explaining why I there may be mixed views on priority.
> > 
> > 
> > 
> > Cheers,
> > 
> > 
> > 
> >   *   Tom
> > 
> > 
> > 
> > From: J3 <j3-bounces at mailman.j3-fortran.org> on behalf of j3 <j3 at mailman.j3-fortran.org>
> > Reply-To: "kargl at uw.edu" <kargl at uw.edu>, j3 <j3 at mailman.j3-fortran.org>
> > Date: Saturday, November 4, 2023 at 11:29 PM
> > To: j3 <j3 at mailman.j3-fortran.org>
> > Cc: "Steven G. Kargl" <kargl at uw.edu>
> > Subject: [EXTERNAL] [BULK] [J3] SCALE(x,i) with complex x
> > 
> > 
> > 
> > CAUTION: This email originated from outside of NASA.  Please take care when clicking links or opening attachments.  Use the "Report Message" button to report suspicious messages to the NASA SOC.
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> > 
> 
> -- 
> Steve
> 

-- 
Steve
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/x-pkcs7-signature
Size: 3950 bytes
Desc: not available
URL: <https://mailman.j3-fortran.org/pipermail/j3/attachments/20231106/a4d9edb2/attachment.bin>