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

[Steven G. Kargl]

On Mon, Nov 06, 2023 at 12:55:20PM +0000, Clune, Thomas L. (GSFC-6101) via J3 wrote:
> 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.

Thanks for the thoughtful reply.   Your guess aligns with what I
also suspect.  Fortran 90 was a massive update to Fortran 77,
and introduced the idea of a model number for the REAL type.  SCALE
was likely motivated by the definition of a model number.

While the code you show is indeed clear to a programmer, it
involves the evaluation of (2**2) (1 multiply), a conversion
of that result to COMPLX(4.,0.), and then the complex multiplication
CMPLX(4.*A%RE-0.*A%IM, 4.*A%IM+0.*A%RE) (4 multiplies, 2 additions).

A Fortran processor might optimize out the multiplications with
0, but this assumes that the result of that multiplication is 
not exceptional.  F2018, 10.1.5.2.4, gives enough wiggle room
to allow this.  An example to consider (2.**2)*(1,Inf) = (NaN,Inf)
while SCALE((1,inf),2) would be (4., Inf).

A more important case might be a point in the 4th quadrant of the
complex (1,sign(0.,-1.)) = (1, -0.) i.e., a point infinitesimally
close to the real axis.  (2**2)*(1., -0.) yields (4., 0.) a point
that is now in the first quadrant.  SCALE((1.,-0.), 2) = (4.,-0.).

-- 
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/fa661b5f/attachment-0001.bin>