Complex arithmetic in Fortran

[Toon Moene]

On 11/13/24 15:40, Richard Biener wrote:
> On Wed, Nov 13, 2024 at 3:21 PM Toon Moene <toon@moene.org> wrote:
>>
>> On 11/13/24 15:12, Richard Biener wrote:
>>
>>> On Wed, Nov 13, 2024 at 3:05 PM Thomas Koenig <tkoenig@netcologne.de> wrote:
>>>>
>>>> Hello world,
>>>>
>>>> J3, the US Fortran standards committee, has passed
>>>> https://j3-fortran.org/doc/year/24/24-179.txt
>>>> which states (with a bit of an overabundance of
>>>> clarity) that, in Fortran, it is possible special-case
>>>> complex multiplication when one of the numbers is known
>>>> to have a zero component, for example when promoting
>>>> a real to complex for complex multiplication.  For
>>>> example, multiplying a complex variable b with a real
>>>> variable a can be done with c%re = b%re * a, c%im = b%im * a,
>>>> without considering NaNs and infinities. Apparently, other
>>>> Fortran compilers do this.
>>>>
>>>> They also stated that ISO/IEC 60559:2020 (aka IEEE 754) does
>>>> not specify complex arithmetic (I wouldn't know, because it is a
>>>> paywalled standard).
>>>>
>>>> How do we want to deal with this? Do we want to implement this
>>>> (it's an obvious speed advantage)?  Should it be the default?
>>>> Do we want to include this in -fcx-fortran-rules?
>>>
>>> The middle-end complex lowering pass does this already, irrespective
>>> of NaNs, same for some degenerate cases with division.
>>
>> Are you sure ?
>>
>> For this code:
>>
>> $ cat complex.f90
>> complex function p(c, r)
>> complex, intent(in) :: c
>> real, intent(in)    :: r
>> p = c * r
>> end
>>
>> I definitely see a difference between
>>
>> $ gfortran -O2 -S complex.f90
>>
>> and
>>
>> $ gfortran -O2 -ffast-math -S complex.f90
> 
> Ah.  This is because of
> 
> static int
> some_nonzerop (tree t)
> {
>    int zerop = false;
> 
>    /* Operations with real or imaginary part of a complex number zero
>       cannot be treated the same as operations with a real or imaginary
>       operand if we care about the signs of zeros in the result.  */
>    if (TREE_CODE (t) == REAL_CST && !flag_signed_zeros)
>      zerop = real_identical (&TREE_REAL_CST (t), &dconst0);
> 
> so the -ffast-math result can be obtained with just -fno-signed-zeros (I assumed
> fortran doesn't have signed zeros?)

Since the Fortran 95 Standard it does (in the current Standard: 7.4.3.2 
Real type):

The real type includes a zero value. Processors that distinguish between 
positive and negative zeros shall treat them as mathematically equivalent
• in all intrinsic relational operations, and
• as actual arguments to intrinsic procedures other than those for which
   it is explicitly specified that negative zero is distinguished.

[Note that "processor" in Fortran standardese means everything 
(combined) from the compiler down to the actual hardware].

So we have to comb through the Standard to see where bullet 2 applies ...

Kind regards,

-- 
Toon Moene - e-mail: toon@moene.org - phone: +31 346 214290
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands