[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" <gcc-bugzilla@gcc.gnu.org>]

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?