[Bug fortran/107753] New: gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[Bug fortran/107753] New: gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[weslley.pereira at ucdenver dot edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

            Bug ID: 107753
           Summary: gfortran returns NaN in complex divisions
                    (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)
           Product: gcc
           Version: 12.2.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: weslley.pereira at ucdenver dot edu
  Target Milestone: ---

If `x=huge(0.0d0)` or `x=2.0d0**(dble(maxexponent(0.0d0))-1)`, the GNU Fortran
12.2.0 returns a NaN for the complex divisions `(x+x*I)/(x+x*I)` and
`(x+x*I)/(x-x*I)`. We verified this after running compiler tests for the new
LAPACK 3.11.0 release. All other divisions with `x=2**m`, for for
`MINEXPONENT-1 <= m < MAXEXPONENT` return the expected results:
`(x+x*I)/(x+x*I)=1` and `(x+x*I)/(x-x*I)=I`.

Related links:
How to reproduce this issue: https://godbolt.org/z/b3WKWodvn
Open issue in LAPACK: https://github.com/Reference-LAPACK/lapack/issues/757
Tests added to LAPACK: https://github.com/Reference-LAPACK/lapack/pull/623

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

anlauf at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
   Last reconfirmed|                            |2022-11-18
             Status|UNCONFIRMED                 |NEW
     Ever confirmed|0                           |1

--- Comment #1 from anlauf at gcc dot gnu.org ---
Confirmed if compiled without optimization (-O0,-Og) on x86_64-pc-linux-gnu.

If I compile with -O, I get:

   (1.79769313486231571E+308,1.79769313486231571E+308)
   (8.98846567431157954E+307,8.98846567431157954E+307)
   (4.49423283715578977E+307,4.49423283715578977E+307)
               (1.0000000000000000,0.0000000000000000)
               (1.0000000000000000,0.0000000000000000)
               (1.0000000000000000,0.0000000000000000)

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org
             Status|NEW                         |WAITING

--- Comment #2 from kargl at gcc dot gnu.org ---
Please either include the program in your original email or attach it the PR.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #3 from anlauf at gcc dot gnu.org ---
I guess the reporter assumes that gcc uses a clever algorithm like Smith's
to handle such extreme cases of complex division.  Not sure if that one is
available by some compilation flag, and I think it would impact performance.

In any case, if the reporter wants to get robust results and in a portable
way, I would advise him to change/fix his algorithm accordingly.  It appears
that a few other compilers behave here like gfortran.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #4 from kargl at gcc dot gnu.org ---
(In reply to anlauf from comment #3)
> I guess the reporter assumes that gcc uses a clever algorithm like Smith's
> to handle such extreme cases of complex division.  Not sure if that one is
> available by some compilation flag, and I think it would impact performance.
> 
> In any case, if the reporter wants to get robust results and in a portable
> way, I would advise him to change/fix his algorithm accordingly.  It appears
> that a few other compilers behave here like gfortran.

It's likely coming from the middle-end where gcc.info has
the option

'-fcx-fortran-rules'
     Complex multiplication and division follow Fortran rules.  Range
     reduction is done as part of complex division, but there is no
     checking whether the result of a complex multiplication or division
     is 'NaN + I*NaN', with an attempt to rescue the situation in that
     case.

Consider

program zdiv
   implicit none
   integer, parameter :: dp = kind(1.d0)
   real(dp) r
   complex(dp) :: y

!   r = huge(1.d0) / 2              ! yields (1,0)
!   r = nearest(huge(1.d0), -1.d0)  ! yields (1,0)
   r = huge(1.d0)                   ! yields (NaN,0)
   y = cmplx(r, r, dp)
   write(*,*) y
   call doit(y)
   write(*,*) y
   contains
      subroutine doit(z)
         complex(dp) z
         z = z / z
      end subroutine doit
end program zdiv

If you compile this with -fdump-tree-all, then one gets
(I've added annotation with <--- marker)


% more z-a.f90.241t.cplxlower0

__attribute__((fn spec (". w ")))
void doit (complex(kind=8) & restrict z)
{
  real(kind=8) D.4265;
  real(kind=8) D.4264;
  complex(kind=8) _1;
  complex(kind=8) _2;
  complex(kind=8) _3;
  real(kind=8) _7;
  real(kind=8) _8;
  real(kind=8) _9;
  real(kind=8) _10;
  real(kind=8) _11;
  real(kind=8) _12;
  logical(kind=1) _13;
  real(kind=8) _14;
  real(kind=8) _15;
  real(kind=8) _16;
  real(kind=8) _17;
  real(kind=8) _18;
  real(kind=8) _19;
  real(kind=8) _20;
  real(kind=8) _21;
  real(kind=8) _22;
  real(kind=8) _23;
  real(kind=8) _24;
  real(kind=8) _25;
  real(kind=8) _26;
  real(kind=8) _27;
  real(kind=8) _28;
  real(kind=8) _29;
  real(kind=8) _30;
  real(kind=8) _31;
  real(kind=8) _32;
  real(kind=8) _33;
  real(kind=8) _34;
  real(kind=8) _35;
  real(kind=8) _36;
  real(kind=8) _37;
  real(kind=8) _38;
  real(kind=8) _39;

  <bb 2> :
  _7 = REALPART_EXPR <*z_5(D)>;
  _8 = IMAGPART_EXPR <*z_5(D)>;
  _1 = COMPLEX_EXPR <_7, _8>;
  _9 = REALPART_EXPR <*z_5(D)>;
  _10 = IMAGPART_EXPR <*z_5(D)>;
  _2 = COMPLEX_EXPR <_9, _10>;
  _11 = ABS_EXPR <_9>;
  _12 = ABS_EXPR <_10>;
  _13 = _11 < _12;
  if (_13 != 0)
    goto <bb 4>; [50.00%]
  else
    goto <bb 5>; [50.00%]

  <bb 4> :
  _14 = _9 / _10;  <--- Should be one
  _15 = _9 * _14;  <--- huge(1.d0)
  _16 = _15 + _10; <--- huge(1.d0) + huge(1.d0) = inf
  _17 = _7 * _14;  <--- huge(1.d0)
  _18 = _17 + _8;  <--- huge(1.d0) + huge(1.d0) = inf
  _19 = _8 * _14;  <--- huge(1.d0)
  _20 = _19 - _7;  <--- huge(1.d0) - huge(1.d0) = 0
  _21 = _18 / _16; <--- inf / inf = NaN
  _22 = _20 / _16; <--- 0 / inf = 0
  _36 = _21;
  _37 = _22;
  goto <bb 3>; [100.00%]

  <bb 5> :
  _23 = _10 / _9;
  _24 = _10 * _23;
  _25 = _24 + _9;
  _26 = _8 * _23;
  _27 = _26 + _7;
  _28 = _7 * _23;
  _29 = _8 - _28;
  _30 = _27 / _25;
  _31 = _29 / _25;
  _38 = _30;
  _39 = _31;

  <bb 3> :
  # _34 = PHI <_36(4), _38(5)>
  # _35 = PHI <_37(4), _39(5)>
  _3 = COMPLEX_EXPR <_34, _35>;
  _32 = _34;
  _33 = _35;
  REALPART_EXPR <*z_5(D)> = _32;
  IMAGPART_EXPR <*z_5(D)> = _33;
  return;

}

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #5 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Fri, Nov 18, 2022 at 10:05:21PM +0000, kargl at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753
> 
> --- Comment #4 from kargl at gcc dot gnu.org ---
> (In reply to anlauf from comment #3)
> > I guess the reporter assumes that gcc uses a clever algorithm like Smith's
> > to handle such extreme cases of complex division.  Not sure if that one is
> > available by some compilation flag, and I think it would impact performance.
> > 
> > In any case, if the reporter wants to get robust results and in a portable
> > way, I would advise him to change/fix his algorithm accordingly.  It appears
> > that a few other compilers behave here like gfortran.
> 
> It's likely coming from the middle-end where gcc.info has
> the option
> 
> '-fcx-fortran-rules'
>      Complex multiplication and division follow Fortran rules.  Range
>      reduction is done as part of complex division, but there is no
>      checking whether the result of a complex multiplication or division
>      is 'NaN + I*NaN', with an attempt to rescue the situation in that
>      case.

Does anyone know what is meant by "Fortran rules"?  F66 does not
have any particular algorithm specified.  I'll look at F77 shortly.

Tracking down what -fcx-fortran-rules does, one finds the
eventually flag_complex_method is set to 1.  The lower of
complex division occurs in gcc/tree-complex.cc (expand_complex_division).
If I use this patch

% git diff gcc/tree-complex.cc | cat
diff --git a/gcc/tree-complex.cc b/gcc/tree-complex.cc
index ea9df6114a1..8051b7a3843 100644
--- a/gcc/tree-complex.cc
+++ b/gcc/tree-complex.cc
@@ -1501,6 +1501,7 @@ expand_complex_division (gimple_stmt_iterator *gsi, tree
type,
          break;

        case 2:
+       case 1:
          if (SCALAR_FLOAT_TYPE_P (inner_type))
            {
              expand_complex_libcall (gsi, type, ar, ai, br, bi, code, true);
@@ -1508,7 +1509,6 @@ expand_complex_division (gimple_stmt_iterator *gsi, tree
type,
            }
          /* FALLTHRU */

-       case 1:
          /* wide ranges of inputs must work for complex divide.  */
          expand_complex_div_wide (gsi, inner_type, ar, ai, br, bi, code);
          break;

to force gfortran through the C language code path, I get

void doit (complex(kind=8) & restrict z)
{
  complex(kind=8) _1;
  complex(kind=8) _2;
  complex(kind=8) _3;
  real(kind=8) _7;
  real(kind=8) _8;
  real(kind=8) _9;
  real(kind=8) _10;
  real(kind=8) _11;
  real(kind=8) _12;

  <bb 2> :
  _7 = REALPART_EXPR <*z_5(D)>;
  _8 = IMAGPART_EXPR <*z_5(D)>;
  _1 = COMPLEX_EXPR <_7, _8>;
  _9 = REALPART_EXPR <*z_5(D)>;
  _10 = IMAGPART_EXPR <*z_5(D)>;
  _2 = COMPLEX_EXPR <_9, _10>;
  _3 = __divdc3 (_7, _8, _9, _10);
  _11 = REALPART_EXPR <_3>;
  _12 = IMAGPART_EXPR <_3>;
  REALPART_EXPR <*z_5(D)> = _11;
  IMAGPART_EXPR <*z_5(D)> = _12;
  return;

}

with the result

%  gfcx -o z -fdump-tree-all a.f90 && ./z
   (1.79769313486231571E+308,1.79769313486231571E+308)
               (1.0000000000000000,0.0000000000000000)

So, is -fcx-fortran-rules a relic of g77 past?

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #6 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Fri, Nov 18, 2022 at 11:24:29PM +0000, sgk at troutmask dot
apl.washington.edu wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753
> 
> --- Comment #5 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
> On Fri, Nov 18, 2022 at 10:05:21PM +0000, kargl at gcc dot gnu.org wrote:
> > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753
> > 
> > --- Comment #4 from kargl at gcc dot gnu.org ---
> > (In reply to anlauf from comment #3)
> > > I guess the reporter assumes that gcc uses a clever algorithm like Smith's
> > > to handle such extreme cases of complex division.  Not sure if that one is
> > > available by some compilation flag, and I think it would impact performance.
> > > 
> > > In any case, if the reporter wants to get robust results and in a portable
> > > way, I would advise him to change/fix his algorithm accordingly.  It appears
> > > that a few other compilers behave here like gfortran.
> > 
> > It's likely coming from the middle-end where gcc.info has
> > the option
> > 
> > '-fcx-fortran-rules'
> >      Complex multiplication and division follow Fortran rules.  Range
> >      reduction is done as part of complex division, but there is no
> >      checking whether the result of a complex multiplication or division
> >      is 'NaN + I*NaN', with an attempt to rescue the situation in that
> >      case.
> 
> Does anyone know what is meant by "Fortran rules"?  F66 does not
> have any particular algorithm specified.  I'll look at F77 shortly.
> 

I add the subroutine

   subroutine ohno
      complex(dp), parameter :: a = cmplx(huge(1.d0),huge(1.d0),dp)
      complex(dp), parameter :: b = a / a
      write(*,*) a
      write(*,*) b
   end subroutine ohno 


% gfortran -o z a.f90 && ./z
   (1.79769313486231571E+308,1.79769313486231571E+308)
                              (NaN,0.0000000000000000)
   (1.79769313486231571E+308,1.79769313486231571E+308)
               (1.0000000000000000,0.0000000000000000)

The last two lines are from ohno.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[weslley.pereira at ucdenver dot edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #7 from Weslley da Silva Pereira <weslley.pereira at ucdenver dot edu> ---
(In reply to anlauf from comment #3)
> I guess the reporter assumes that gcc uses a clever algorithm like Smith's
> to handle such extreme cases of complex division.  Not sure if that one is
> available by some compilation flag, and I think it would impact performance.
> 
> In any case, if the reporter wants to get robust results and in a portable
> way, I would advise him to change/fix his algorithm accordingly.  It appears
> that a few other compilers behave here like gfortran.

Thanks for the suggestion of changing the algorithm that needs such a division.
What are the ranges for nominator and denominator where one should rely on the
intrinsic complex division? Maybe this is the good question to ask. Then we can
build algorithms that attend such requisites. LAPACK has cladiv and zladiv,
which are routines for complex division that avoids unnecessary under and
overflow. They are used in many parts of the code. This implementation is for
sure less efficient than the intrinsic complex division, but we rely on it
because of robustness.

More data for the discussion:
1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization
flags `-O` but still reproduce the wrong result for complex divisions with huge
numbers. See
https://github.com/Reference-LAPACK/lapack/issues/575#issuecomment-910616816
that used the code from
https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/test_zcomplexdiv.f.
This is the test currently in LAPACK 3.11.0.
2. I have just reproduced what was reported in
https://github.com/Reference-LAPACK/lapack/issues/575#issuecomment-910616816 in
my Ubuntu 20.04.5 LTS, using GNU Fortran 9.4.0.
3. I noticed that the optimization flag is unable to target divisions like
`x/x` depending on where they are inside a program.
4. My Ubuntu 20.04.5 LTS with compiler ifort 2021.7.1 computes the complex
division `x/x` accurately even for the case of huge numbers. Scenarios tested:
   - I tested the program in
https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/test_zcomplexdiv.f
and the one in https://godbolt.org/z/b3WKWodvn.
   - I tested ifort with flags -fp-model precise and -fp-model fast. The latter
enables more aggressive optimizations on floating-point data.
   - I tested compilation with optimization flags -O0, -O, -O1, -O2, -O3. 

Here is the implementation of the complex division in LAPACK if it somehow
helps the discussion:
https://netlib.org/lapack/explore-html/d8/d9b/group__double_o_t_h_e_rauxiliary_gad1c0279ec29e8ac222f1e319f4144fcb.html

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[weslley.pereira at ucdenver dot edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #8 from Weslley da Silva Pereira <weslley.pereira at ucdenver dot edu> ---
Created attachment 53927
  --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=53927&action=edit
Test case with many examples of complex division

Test code used in LAPACK 3.11.0.
Code extracted from
https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/test_zcomplexdiv.f

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #9 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Fri, Nov 18, 2022 at 11:24:29PM +0000, sgk at troutmask dot
apl.washington.edu wrote:
> 
> Does anyone know what is meant by "Fortran rules"?  F66 does not
> have any particular algorithm specified.  I'll look at F77 shortly.
> 

Well, I hunted down the origins of -fcx-fortran-rules.

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=29549

So, it appears to be an optimization, where Smith's algorithm
will fail for extreme values of the real and imaginary parts
of the complex number.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Priority|P3                          |P4

--- Comment #10 from kargl at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #9)
> On Fri, Nov 18, 2022 at 11:24:29PM +0000, sgk at troutmask dot
> apl.washington.edu wrote:
> > 
> > Does anyone know what is meant by "Fortran rules"?  F66 does not
> > have any particular algorithm specified.  I'll look at F77 shortly.
> > 
> 
> Well, I hunted down the origins of -fcx-fortran-rules.
> 
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=29549
> 
> So, it appears to be an optimization, where Smith's algorithm
> will fail for extreme values of the real and imaginary parts
> of the complex number.

So, I wrong a dirty little program to time complex division.


program foo

   use timerm, only : rdtsc

   implicit none

   integer, parameter :: n = 1024*1024, dp = kind(1.d0)
   real(dp) re(n), im(n)
   complex(dp) x(n)
   integer i
   integer(8) t0, t

   t = 0
   do i = 1, 10
      call random_number(re)
      call random_number(im)
      x = cmplx(4 * re, 10 * im,8)
      t0 = rdtsc()
      x = x / x
      t = t + (rdtsc() - t0)
   end do
   print '(G0,1X,G0)', x(1)
   print *, real(t,10) / 10 / n
end program foo


Compiled with gfortran with its current method of doing division
(i.e., -fcx-fortran-rules), I see roughly 44.5 clock ticks per
division.  If run with a patched gfortran that uses the method that
the C compiler uses, I get about 62 ticks per division.  So, using
the stricter method impacts performance.

I'll note that gfortran unilaterally enforces -fcx-fortran-rules,
i.e., -fno-cx-fortran-rules has no effect.  Perhaps, gfortran 
could be given a new -fcx-division=XXX option, where XXX is one of

   naive --> what -ffast-math does   (flags_complex_method = 0)
   smith --> what -fcx-fortran-rules (flags_complex_method = 1)
   strict -> default C behavior      (flags_complex_method = 2)

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #11 from anlauf at gcc dot gnu.org ---
(In reply to Weslley da Silva Pereira from comment #7)
> More data for the discussion:
> 1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization
> flags `-O` but still reproduce the wrong result for complex divisions with
> huge numbers. See

It is possible that gfortran's dependence on optimization level depends
on the version.  If one wants to test run-time behavior and avoid
compile-time simplification, it may be helpful to add:

  volatile :: x, y, z

I then get consistent results for -O0 / -O1.

> 4. My Ubuntu 20.04.5 LTS with compiler ifort 2021.7.1 computes the complex
> division `x/x` accurately even for the case of huge numbers. Scenarios
> tested:
>    - I tested the program in
> https://github.com/Reference-LAPACK/lapack/blob/master/INSTALL/
> test_zcomplexdiv.f and the one in https://godbolt.org/z/b3WKWodvn.
>    - I tested ifort with flags -fp-model precise and -fp-model fast. The
> latter enables more aggressive optimizations on floating-point data.
>    - I tested compilation with optimization flags -O0, -O, -O1, -O2, -O3. 

Intel might be fine, but at least some current llvm-based compilers (Nvidia,
AMD flang) show more or less similar behavior to gfortran.

E.g. nvfortran 22.11:

 (1.7976931348623157E+308,1.7976931348623157E+308)
 (8.9884656743115795E+307,8.9884656743115795E+307)
 (4.4942328371557898E+307,4.4942328371557898E+307)
 (NaN,0.000000000000000)
 (NaN,0.000000000000000)
 (1.000000000000000,0.000000000000000)

As a sidenote: we are really discussing borderline cases here, valid
but only rarely occuring in normal code execution.  If I replace

  x = cmplx( huge(0.0d0), huge(0.0d0), dp )
  y = cmplx( b**(E-1), b**(E-1), dp )

by

  x = cmplx( nearest(huge(0.0d0),-1.d0), nearest(huge(0.0d0),-1.d0), dp )
  y = cmplx( nearest(b**(E-1),   -1.d0), nearest(b**(E-1),   -1.d0), dp )

then I get

               (1.0000000000000000,0.0000000000000000)
               (1.0000000000000000,0.0000000000000000)

instead of NaN.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #12 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Sat, Nov 19, 2022 at 08:14:01PM +0000, anlauf at gcc dot gnu.org wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753
> 
> --- Comment #11 from anlauf at gcc dot gnu.org ---
> (In reply to Weslley da Silva Pereira from comment #7)
> > More data for the discussion:
> > 1. In a Ubuntu 18.04.5 LTS, using GNU Fortran 7.5.0, I tested optimization
> > flags `-O` but still reproduce the wrong result for complex divisions with
> > huge numbers. See
> 
> It is possible that gfortran's dependence on optimization level depends
> on the version.  If one wants to test run-time behavior and avoid
> compile-time simplification, it may be helpful to add:
> 
>   volatile :: x, y, z
> 
> I then get consistent results for -O0 / -O1.
> 

The optimization level is irrelevant.  gfortran unilaterally
uses -fcx-fortran-rules, and there is no way to disable this
option to user the slower, but stricter, evaluation.  One
will always get complex division computed by

a+ib   a + b(d/c)     b - a(d/c) 
---- = ---------- + i ------------  |c| > |d|
c+id   c + d(d/c)     c + d(d/c)

and similar for |d| > |c|.

There are a few problems with this. d/c can trigger an invalid underflow
exception.  If d == c, you then have numerators of a + b and b - a, you
can get a invalid overflow for a = huge() and b > 1e291_8.

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #13 from anlauf at gcc dot gnu.org ---
(In reply to Steve Kargl from comment #12)
> The optimization level is irrelevant.  gfortran unilaterally
> uses -fcx-fortran-rules, and there is no way to disable this
> option to user the slower, but stricter, evaluation.  One
> will always get complex division computed by
> 
> a+ib   a + b(d/c)     b - a(d/c) 
> ---- = ---------- + i ------------  |c| > |d|
> c+id   c + d(d/c)     c + d(d/c)
> 
> and similar for |d| > |c|.
> 
> There are a few problems with this. d/c can trigger an invalid underflow
> exception.  If d == c, you then have numerators of a + b and b - a, you
> can get a invalid overflow for a = huge() and b > 1e291_8.

I am wondering how slow an algorithm would be that scales numerator
and denominator by respective factors that are powers of 2, e.g.

e_num = 2. ** -max (exponent (a), exponent (b))
e_den = 2. ** -max (exponent (c), exponent (d))

The modulus of scaled values would be <= 1, even for any of a,... being huge().
Of course this does not address underflows that could occur during scaling,
or denormalized numbers, which are numerically irrelevant for the result.

Is there anything else wrong with this approach?

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #14 from kargl at gcc dot gnu.org ---
(In reply to anlauf from comment #13)
> (In reply to Steve Kargl from comment #12)
> > The optimization level is irrelevant.  gfortran unilaterally
> > uses -fcx-fortran-rules, and there is no way to disable this
> > option to user the slower, but stricter, evaluation.  One
> > will always get complex division computed by
> > 
> > a+ib   a + b(d/c)     b - a(d/c) 
> > ---- = ---------- + i ------------  |c| > |d|
> > c+id   c + d(d/c)     c + d(d/c)
> > 
> > and similar for |d| > |c|.
> > 
> > There are a few problems with this. d/c can trigger an invalid underflow
> > exception.  If d == c, you then have numerators of a + b and b - a, you
> > can get a invalid overflow for a = huge() and b > 1e291_8.
> 
> I am wondering how slow an algorithm would be that scales numerator
> and denominator by respective factors that are powers of 2, e.g.
> 
> e_num = 2. ** -max (exponent (a), exponent (b))
> e_den = 2. ** -max (exponent (c), exponent (d))
> 
> The modulus of scaled values would be <= 1, even for any of a,... being
> huge().
> Of course this does not address underflows that could occur during scaling,
> or denormalized numbers, which are numerically irrelevant for the result.
> 
> Is there anything else wrong with this approach?

Comment #10 contains a simple timing measurement in from my Intel Core2 Duo
based system.  gfortran with its current method (ie., -fcx-fortran-rules) takes
44.5 clock ticks for a complex division.  If I sidestep the option and force it
to use the C language method of evaluation, it takes 62 clock ticks.  I haven't
looked at what algorithm C uses, but I suspect its along the lines you suggest.
 The question is likely do we break backwards compatibility and remove
-fcx-fortran-rules or change when/how -fcx-fortran-rules applies (e.g., add it
to -ffast-math?)

[Bug fortran/107753] gfortran returns NaN in complex divisions (x+x*I)/(x+x*I) and (x+x*I)/(x-x*I)

[anlauf at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=107753

--- Comment #15 from anlauf at gcc dot gnu.org ---
(In reply to kargl from comment #14)
> Comment #10 contains a simple timing measurement in from my Intel Core2 Duo
> based system.  gfortran with its current method (ie., -fcx-fortran-rules)
> takes 44.5 clock ticks for a complex division.
[...]
> ticks.  I haven't looked at what algorithm C uses, but I suspect its along
> the lines you suggest.

I may have misread the library implementation of complex division (libgcc2),
but it likely does really much more that we need for Fortran.  And some
of the overhead seems to come from doing many more divisions in the
library version than in the current gfortran variant.

> The question is likely do we break backwards
> compatibility and remove -fcx-fortran-rules or change when/how
> -fcx-fortran-rules applies (e.g., add it to -ffast-math?)

I would certainly not be happy breaking backwards compatibility,
but is there a precise definition of what the "cx-fortran-rules" are?

I had the silent hope that manipulating exponents of IEEE floats could be a
cheap operation, allowing for a cheaper solution than the expensive method
that gets really borderline cases better, although an unreasonable cost.